- Timestamp:
- Apr 9, 2010, 8:49:12 PM (14 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
anuga_core/source/anuga/abstract_2d_finite_volumes/util.py
r7675 r7679 15 15 16 16 from anuga.utilities.numerical_tools import ensure_numeric 17 from Scientific.IO.NetCDF import NetCDFFile 18 19 from anuga.geospatial_data.geospatial_data import ensure_absolute 17 20 18 from math import sqrt, atan, degrees 21 19 … … 25 23 from anuga.utilities.system_tools import store_version_info 26 24 27 from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a28 25 import anuga.utilities.log as log 29 26 … … 31 28 32 29 33 ##34 # @brief Read time history of data from NetCDF file, return callable object.35 # @param filename Name of .sww or .tms file.36 # @param domain Associated domain object.37 # @param quantities Name of quantity to be interpolated or a list of names.38 # @param interpolation_points List of absolute UTM coordinates for points39 # (N x 2) or geospatial object or40 # points file name at which values are sought.41 # @param time_thinning42 # @param verbose True if this function is to be verbose.43 # @param use_cache True means that caching of intermediate result is attempted.44 # @param boundary_polygon45 # @param output_centroids if True, data for the centroid of the triangle will be output46 # @return A callable object.47 30 def file_function(filename, 48 31 domain=None, … … 55 38 boundary_polygon=None, 56 39 output_centroids=False): 57 """Read time history of spatial data from NetCDF file and return 58 a callable object. 59 60 Input variables: 61 62 filename - Name of sww, tms or sts file 63 64 If the file has extension 'sww' then it is assumed to be spatio-temporal 65 or temporal and the callable object will have the form f(t,x,y) or f(t) 66 depending on whether the file contains spatial data 67 68 If the file has extension 'tms' then it is assumed to be temporal only 69 and the callable object will have the form f(t) 70 71 Either form will return interpolated values based on the input file 72 using the underlying interpolation_function. 73 74 domain - Associated domain object 75 If domain is specified, model time (domain.starttime) 76 will be checked and possibly modified. 77 78 All times are assumed to be in UTC 79 80 All spatial information is assumed to be in absolute UTM coordinates. 81 82 quantities - the name of the quantity to be interpolated or a 83 list of quantity names. The resulting function will return 84 a tuple of values - one for each quantity 85 If quantities are None, the default quantities are 86 ['stage', 'xmomentum', 'ymomentum'] 87 88 89 interpolation_points - list of absolute UTM coordinates for points (N x 2) 90 or geospatial object or points file name at which values are sought 91 92 time_thinning - 93 94 verbose - 95 96 use_cache: True means that caching of intermediate result of 97 Interpolation_function is attempted 98 99 boundary_polygon - 100 101 102 See Interpolation function in anuga.fit_interpolate.interpolation for 103 further documentation 104 """ 105 106 # FIXME (OLE): Should check origin of domain against that of file 107 # In fact, this is where origin should be converted to that of domain 108 # Also, check that file covers domain fully. 109 110 # Take into account: 111 # - domain's georef 112 # - sww file's georef 113 # - interpolation points as absolute UTM coordinates 114 115 if quantities is None: 116 if verbose: 117 msg = 'Quantities specified in file_function are None,' 118 msg += ' so I will use stage, xmomentum, and ymomentum in that order' 119 log.critical(msg) 120 quantities = ['stage', 'xmomentum', 'ymomentum'] 121 122 # Use domain's startime if available 123 if domain is not None: 124 domain_starttime = domain.get_starttime() 125 else: 126 domain_starttime = None 127 128 # Build arguments and keyword arguments for use with caching or apply. 129 args = (filename,) 130 131 # FIXME (Ole): Caching this function will not work well 132 # if domain is passed in as instances change hash code. 133 # Instead we pass in those attributes that are needed (and return them 134 # if modified) 135 kwargs = {'quantities': quantities, 136 'interpolation_points': interpolation_points, 137 'domain_starttime': domain_starttime, 138 'time_thinning': time_thinning, 139 'time_limit': time_limit, 140 'verbose': verbose, 141 'boundary_polygon': boundary_polygon, 142 'output_centroids': output_centroids} 143 144 # Call underlying engine with or without caching 145 if use_cache is True: 146 try: 147 from caching import cache 148 except: 149 msg = 'Caching was requested, but caching module'+\ 150 'could not be imported' 151 raise msg 152 153 f, starttime = cache(_file_function, 154 args, kwargs, 155 dependencies=[filename], 156 compression=False, 157 verbose=verbose) 158 else: 159 f, starttime = apply(_file_function, 160 args, kwargs) 161 162 #FIXME (Ole): Pass cache arguments, such as compression, in some sort of 163 #structure 164 165 f.starttime = starttime 166 f.filename = filename 167 168 if domain is not None: 169 #Update domain.startime if it is *earlier* than starttime from file 170 if starttime > domain.starttime: 171 msg = 'WARNING: Start time as specified in domain (%f)'\ 172 %domain.starttime 173 msg += ' is earlier than the starttime of file %s (%f).'\ 174 %(filename, starttime) 175 msg += ' Modifying domain starttime accordingly.' 176 177 if verbose: log.critical(msg) 178 179 domain.set_starttime(starttime) #Modifying model time 180 181 if verbose: log.critical('Domain starttime is now set to %f' 182 % domain.starttime) 183 return f 184 185 186 ## 187 # @brief ?? 188 # @param filename Name of .sww or .tms file. 189 # @param domain Associated domain object. 190 # @param quantities Name of quantity to be interpolated or a list of names. 191 # @param interpolation_points List of absolute UTM coordinates for points 192 # (N x 2) or geospatial object or 193 # points file name at which values are sought. 194 # @param time_thinning 195 # @param verbose True if this function is to be verbose. 196 # @param use_cache True means that caching of intermediate result is attempted. 197 # @param boundary_polygon 198 def _file_function(filename, 199 quantities=None, 200 interpolation_points=None, 201 domain_starttime=None, 202 time_thinning=1, 203 time_limit=None, 204 verbose=False, 205 boundary_polygon=None, 206 output_centroids=False): 207 """Internal function 208 209 See file_function for documentatiton 210 """ 211 212 assert type(filename) == type(''),\ 213 'First argument to File_function must be a string' 214 215 try: 216 fid = open(filename) 217 except Exception, e: 218 msg = 'File "%s" could not be opened: Error="%s"' % (filename, e) 219 raise msg 220 221 # read first line of file, guess file type 222 line = fid.readline() 223 fid.close() 224 225 if line[:3] == 'CDF': 226 return get_netcdf_file_function(filename, 227 quantities, 228 interpolation_points, 229 domain_starttime, 230 time_thinning=time_thinning, 231 time_limit=time_limit, 232 verbose=verbose, 233 boundary_polygon=boundary_polygon, 234 output_centroids=output_centroids) 235 else: 236 # FIXME (Ole): Could add csv file here to address Ted Rigby's 237 # suggestion about reading hydrographs. 238 # This may also deal with the gist of ticket:289 239 raise 'Must be a NetCDF File' 240 241 242 ## 243 # @brief ?? 244 # @param filename Name of .sww or .tms file. 245 # @param quantity_names Name of quantity to be interpolated or a list of names. 246 # @param interpolation_points List of absolute UTM coordinates for points 247 # (N x 2) or geospatial object or 248 # points file name at which values are sought. 249 # @param domain_starttime Start time from domain object. 250 # @param time_thinning ?? 251 # @param verbose True if this function is to be verbose. 252 # @param boundary_polygon ?? 253 # @return A callable object. 254 def get_netcdf_file_function(filename, 255 quantity_names=None, 256 interpolation_points=None, 257 domain_starttime=None, 258 time_thinning=1, 259 time_limit=None, 260 verbose=False, 261 boundary_polygon=None, 262 output_centroids=False): 263 """Read time history of spatial data from NetCDF sww file and 264 return a callable object f(t,x,y) 265 which will return interpolated values based on the input file. 266 267 Model time (domain_starttime) 268 will be checked, possibly modified and returned 269 270 All times are assumed to be in UTC 271 272 See Interpolation function for further documetation 273 """ 274 275 # FIXME: Check that model origin is the same as file's origin 276 # (both in UTM coordinates) 277 # If not - modify those from file to match domain 278 # (origin should be passed in) 279 # Take this code from e.g. dem2pts in data_manager.py 280 # FIXME: Use geo_reference to read and write xllcorner... 281 282 import time, calendar, types 283 from anuga.config import time_format 284 285 # Open NetCDF file 286 if verbose: log.critical('Reading %s' % filename) 287 288 fid = NetCDFFile(filename, netcdf_mode_r) 289 290 if type(quantity_names) == types.StringType: 291 quantity_names = [quantity_names] 292 293 if quantity_names is None or len(quantity_names) < 1: 294 msg = 'No quantities are specified in file_function' 295 raise Exception, msg 296 297 if interpolation_points is not None: 298 interpolation_points = ensure_absolute(interpolation_points) 299 msg = 'Points must by N x 2. I got %d' % interpolation_points.shape[1] 300 assert interpolation_points.shape[1] == 2, msg 301 302 # Now assert that requested quantitites (and the independent ones) 303 # are present in file 304 missing = [] 305 for quantity in ['time'] + quantity_names: 306 if not fid.variables.has_key(quantity): 307 missing.append(quantity) 308 309 if len(missing) > 0: 310 msg = 'Quantities %s could not be found in file %s'\ 311 % (str(missing), filename) 312 fid.close() 313 raise Exception, msg 314 315 # Decide whether this data has a spatial dimension 316 spatial = True 317 for quantity in ['x', 'y']: 318 if not fid.variables.has_key(quantity): 319 spatial = False 320 321 if filename[-3:] == 'tms' and spatial is True: 322 msg = 'Files of type tms must not contain spatial information' 323 raise msg 324 325 if filename[-3:] == 'sww' and spatial is False: 326 msg = 'Files of type sww must contain spatial information' 327 raise msg 328 329 if filename[-3:] == 'sts' and spatial is False: 330 #What if mux file only contains one point 331 msg = 'Files of type sts must contain spatial information' 332 raise msg 333 334 if filename[-3:] == 'sts' and boundary_polygon is None: 335 #What if mux file only contains one point 336 msg = 'Files of type sts require boundary polygon' 337 raise msg 338 339 # Get first timestep 340 try: 341 starttime = fid.starttime[0] 342 except ValueError: 343 msg = 'Could not read starttime from file %s' % filename 344 raise msg 345 346 # Get variables 347 # if verbose: log.critical('Get variables' ) 348 time = fid.variables['time'][:] 349 # FIXME(Ole): Is time monotoneous? 350 351 # Apply time limit if requested 352 upper_time_index = len(time) 353 msg = 'Time vector obtained from file %s has length 0' % filename 354 assert upper_time_index > 0, msg 355 356 if time_limit is not None: 357 # Adjust given time limit to given start time 358 time_limit = time_limit - starttime 359 360 361 # Find limit point 362 for i, t in enumerate(time): 363 if t > time_limit: 364 upper_time_index = i 365 break 366 367 msg = 'Time vector is zero. Requested time limit is %f' % time_limit 368 assert upper_time_index > 0, msg 369 370 if time_limit < time[-1] and verbose is True: 371 log.critical('Limited time vector from %.2fs to %.2fs' 372 % (time[-1], time_limit)) 373 374 time = time[:upper_time_index] 375 376 377 378 379 # Get time independent stuff 380 if spatial: 381 # Get origin 382 xllcorner = fid.xllcorner[0] 383 yllcorner = fid.yllcorner[0] 384 zone = fid.zone[0] 385 386 x = fid.variables['x'][:] 387 y = fid.variables['y'][:] 388 if filename.endswith('sww'): 389 triangles = fid.variables['volumes'][:] 390 391 x = num.reshape(x, (len(x),1)) 392 y = num.reshape(y, (len(y),1)) 393 vertex_coordinates = num.concatenate((x,y), axis=1) #m x 2 array 394 395 if boundary_polygon is not None: 396 # Remove sts points that do not lie on boundary 397 # FIXME(Ole): Why don't we just remove such points from the list of points and associated data? 398 # I am actually convinced we can get rid of neighbour_gauge_id altogether as the sts file is produced using the ordering file. 399 # All sts points are therefore always present in the boundary. In fact, they *define* parts of the boundary. 400 boundary_polygon=ensure_numeric(boundary_polygon) 401 boundary_polygon[:,0] -= xllcorner 402 boundary_polygon[:,1] -= yllcorner 403 temp=[] 404 boundary_id=[] 405 gauge_id=[] 406 for i in range(len(boundary_polygon)): 407 for j in range(len(x)): 408 if num.allclose(vertex_coordinates[j], 409 boundary_polygon[i], 1e-4): 410 #FIXME: 411 #currently gauges lat and long is stored as float and 412 #then cast to double. This cuases slight repositioning 413 #of vertex_coordinates. 414 temp.append(boundary_polygon[i]) 415 gauge_id.append(j) 416 boundary_id.append(i) 417 break 418 gauge_neighbour_id=[] 419 for i in range(len(boundary_id)-1): 420 if boundary_id[i]+1==boundary_id[i+1]: 421 gauge_neighbour_id.append(i+1) 422 else: 423 gauge_neighbour_id.append(-1) 424 if boundary_id[len(boundary_id)-1]==len(boundary_polygon)-1 \ 425 and boundary_id[0]==0: 426 gauge_neighbour_id.append(0) 427 else: 428 gauge_neighbour_id.append(-1) 429 gauge_neighbour_id=ensure_numeric(gauge_neighbour_id) 430 431 if len(num.compress(gauge_neighbour_id>=0,gauge_neighbour_id)) \ 432 != len(temp)-1: 433 msg='incorrect number of segments' 434 raise msg 435 vertex_coordinates=ensure_numeric(temp) 436 if len(vertex_coordinates)==0: 437 msg = 'None of the sts gauges fall on the boundary' 438 raise msg 439 else: 440 gauge_neighbour_id=None 441 442 if interpolation_points is not None: 443 # Adjust for georef 444 interpolation_points[:,0] -= xllcorner 445 interpolation_points[:,1] -= yllcorner 446 else: 447 gauge_neighbour_id=None 448 449 if domain_starttime is not None: 450 # If domain_startime is *later* than starttime, 451 # move time back - relative to domain's time 452 if domain_starttime > starttime: 453 time = time - domain_starttime + starttime 454 455 # FIXME Use method in geo to reconcile 456 # if spatial: 457 # assert domain.geo_reference.xllcorner == xllcorner 458 # assert domain.geo_reference.yllcorner == yllcorner 459 # assert domain.geo_reference.zone == zone 460 461 if verbose: 462 log.critical('File_function data obtained from: %s' % filename) 463 log.critical(' References:') 464 if spatial: 465 log.critical(' Lower left corner: [%f, %f]' 466 % (xllcorner, yllcorner)) 467 log.critical(' Start time: %f' % starttime) 468 469 470 # Produce values for desired data points at 471 # each timestep for each quantity 472 quantities = {} 473 for i, name in enumerate(quantity_names): 474 quantities[name] = fid.variables[name][:] 475 if boundary_polygon is not None: 476 #removes sts points that do not lie on boundary 477 quantities[name] = num.take(quantities[name], gauge_id, axis=1) 478 479 # Close sww, tms or sts netcdf file 480 fid.close() 481 482 from anuga.fit_interpolate.interpolate import Interpolation_function 483 484 if not spatial: 485 vertex_coordinates = triangles = interpolation_points = None 486 if filename[-3:] == 'sts':#added 487 triangles = None 488 #vertex coordinates is position of urs gauges 489 490 if verbose: 491 log.critical('Calling interpolation function') 492 493 # Return Interpolation_function instance as well as 494 # starttime for use to possible modify that of domain 495 return (Interpolation_function(time, 496 quantities, 497 quantity_names, 498 vertex_coordinates, 499 triangles, 500 interpolation_points, 501 time_thinning=time_thinning, 502 verbose=verbose, 503 gauge_neighbour_id=gauge_neighbour_id, 504 output_centroids=output_centroids), 505 starttime) 506 507 # NOTE (Ole): Caching Interpolation function is too slow as 508 # the very long parameters need to be hashed. 509 40 from file_function import file_function as file_function_new 41 return file_function_new(filename, domain, quantities, interpolation_points, 42 time_thinning, time_limit, verbose, use_cache, 43 boundary_polygon, output_centroids) 510 44 511 45 ## … … 645 179 return NT.anglediff(v0, v1) 646 180 647 648 # @note TEMP649 def mean(x):650 """Temporary Interface to new location"""651 652 import anuga.utilities.numerical_tools as NT653 654 msg = 'mean has moved from util.py. '655 msg += 'Please use "from anuga.utilities.numerical_tools import mean"'656 warn(msg, DeprecationWarning)657 658 return NT.mean(x)659 660 661 181 # @note TEMP 662 182 def point_on_line(*args, **kwargs): … … 744 264 745 265 746 ##747 # @brief Read a .sww file and plot the time series.748 # @note This function is deprecated - use gauge.sww2timeseries instead.749 #750 def sww2timeseries(swwfiles,751 gauge_filename,752 production_dirs,753 report=None,754 reportname=None,755 plot_quantity=None,756 generate_fig=False,757 surface=None,758 time_min=None,759 time_max=None,760 time_thinning=1,761 time_unit=None,762 title_on=None,763 use_cache=False,764 verbose=False):765 return gauge.sww2timeseries(swwfiles, gauge_filename, production_dirs, report, reportname, \766 plot_quantity, generate_fig, surface, time_min, time_max, \767 time_thinning, time_unit, title_on, use_cache, verbose)768 769 770 266 771 267 ## … … 774 270 # @return A (gauges, gaugelocation, elev) tuple. 775 271 def get_gauges_from_file(filename): 776 return gauge .get_from_file(filename)272 return gauge_get_from_file(filename) 777 273 778 274 … … 2078 1574 This is really returning speed, not velocity. 2079 1575 """ 2080 from gauge import sww2csv 1576 from gauge import sww2csv_gauges as sww2csv 2081 1577 2082 sww2csv(sww_file, gauge_file, out_name, quantities, verbose, use_cache)1578 return sww2csv(sww_file, gauge_file, out_name, quantities, verbose, use_cache) 2083 1579 1580 def sww2timeseries(swwfiles, 1581 gauge_filename, 1582 production_dirs, 1583 report=None, 1584 reportname=None, 1585 plot_quantity=None, 1586 generate_fig=False, 1587 surface=None, 1588 time_min=None, 1589 time_max=None, 1590 time_thinning=1, 1591 time_unit=None, 1592 title_on=None, 1593 use_cache=False, 1594 verbose=False, 1595 output_centroids=False): 1596 from gauge import sww2timeseries as sww2timeseries_new 1597 return sww2timeseries_new(swwfiles, 1598 gauge_filename, 1599 production_dirs, 1600 report, 1601 reportname, 1602 plot_quantity, 1603 generate_fig, 1604 surface, 1605 time_min, 1606 time_max, 1607 time_thinning, 1608 time_unit, 1609 title_on, 1610 use_cache, 1611 verbose, 1612 output_centroids) 2084 1613 2085 1614 ##
Note: See TracChangeset
for help on using the changeset viewer.