[2792] | 1 | """datamanager.py - input output for AnuGA |
---|
| 2 | |
---|
| 3 | |
---|
| 4 | This module takes care of reading and writing datafiles such as topograhies, |
---|
| 5 | model output, etc |
---|
| 6 | |
---|
| 7 | |
---|
| 8 | Formats used within AnuGA: |
---|
| 9 | |
---|
| 10 | .sww: Netcdf format for storing model output f(t,x,y) |
---|
| 11 | .tms: Netcdf format for storing time series f(t) |
---|
| 12 | |
---|
| 13 | .xya: ASCII format for storing arbitrary points and associated attributes |
---|
| 14 | .pts: NetCDF format for storing arbitrary points and associated attributes |
---|
| 15 | |
---|
| 16 | .asc: ASCII format of regular DEMs as output from ArcView |
---|
| 17 | .prj: Associated ArcView file giving more meta data for asc format |
---|
| 18 | .ers: ERMapper header format of regular DEMs for ArcView |
---|
| 19 | |
---|
| 20 | .dem: NetCDF representation of regular DEM data |
---|
| 21 | |
---|
| 22 | .tsh: ASCII format for storing meshes and associated boundary and region info |
---|
| 23 | .msh: NetCDF format for storing meshes and associated boundary and region info |
---|
| 24 | |
---|
| 25 | .nc: Native ferret NetCDF format |
---|
| 26 | .geo: Houdinis ascii geometry format (?) |
---|
| 27 | |
---|
| 28 | |
---|
| 29 | A typical dataflow can be described as follows |
---|
| 30 | |
---|
| 31 | Manually created files: |
---|
| 32 | ASC, PRJ: Digital elevation models (gridded) |
---|
| 33 | TSH: Triangular meshes (e.g. created from pmesh) |
---|
| 34 | NC Model outputs for use as boundary conditions (e.g from MOST) |
---|
| 35 | |
---|
| 36 | |
---|
| 37 | AUTOMATICALLY CREATED FILES: |
---|
| 38 | |
---|
| 39 | ASC, PRJ -> DEM -> PTS: Conversion of DEM's to native pts file |
---|
| 40 | |
---|
| 41 | NC -> SWW: Conversion of MOST bundary files to boundary sww |
---|
| 42 | |
---|
| 43 | PTS + TSH -> TSH with elevation: Least squares fit |
---|
| 44 | |
---|
| 45 | TSH -> SWW: Conversion of TSH to sww viewable using Swollen |
---|
| 46 | |
---|
| 47 | TSH + Boundary SWW -> SWW: Simluation using pyvolution |
---|
| 48 | |
---|
| 49 | """ |
---|
| 50 | import os |
---|
| 51 | |
---|
| 52 | from Numeric import concatenate, array, Float, Int, Int32, resize, sometrue, \ |
---|
| 53 | searchsorted, zeros, allclose, around |
---|
| 54 | |
---|
| 55 | |
---|
| 56 | from coordinate_transforms.geo_reference import Geo_reference |
---|
| 57 | |
---|
| 58 | def make_filename(s): |
---|
| 59 | """Transform argument string into a suitable filename |
---|
| 60 | """ |
---|
| 61 | |
---|
| 62 | s = s.strip() |
---|
| 63 | s = s.replace(' ', '_') |
---|
| 64 | s = s.replace('(', '') |
---|
| 65 | s = s.replace(')', '') |
---|
| 66 | s = s.replace('__', '_') |
---|
| 67 | |
---|
| 68 | return s |
---|
| 69 | |
---|
| 70 | |
---|
| 71 | def check_dir(path, verbose=None): |
---|
| 72 | """Check that specified path exists. |
---|
| 73 | If path does not exist it will be created if possible |
---|
| 74 | |
---|
| 75 | USAGE: |
---|
| 76 | checkdir(path, verbose): |
---|
| 77 | |
---|
| 78 | ARGUMENTS: |
---|
| 79 | path -- Directory |
---|
| 80 | verbose -- Flag verbose output (default: None) |
---|
| 81 | |
---|
| 82 | RETURN VALUE: |
---|
| 83 | Verified path including trailing separator |
---|
| 84 | |
---|
| 85 | """ |
---|
| 86 | |
---|
| 87 | import os, sys |
---|
| 88 | import os.path |
---|
| 89 | |
---|
| 90 | if sys.platform in ['nt', 'dos', 'win32', 'what else?']: |
---|
| 91 | unix = 0 |
---|
| 92 | else: |
---|
| 93 | unix = 1 |
---|
| 94 | |
---|
| 95 | |
---|
| 96 | if path[-1] != os.sep: |
---|
| 97 | path = path + os.sep # Add separator for directories |
---|
| 98 | |
---|
| 99 | path = os.path.expanduser(path) # Expand ~ or ~user in pathname |
---|
| 100 | if not (os.access(path,os.R_OK and os.W_OK) or path == ''): |
---|
| 101 | try: |
---|
| 102 | exitcode=os.mkdir(path) |
---|
| 103 | |
---|
| 104 | # Change access rights if possible |
---|
| 105 | # |
---|
| 106 | if unix: |
---|
| 107 | exitcode=os.system('chmod 775 '+path) |
---|
| 108 | else: |
---|
| 109 | pass # FIXME: What about acces rights under Windows? |
---|
| 110 | |
---|
| 111 | if verbose: print 'MESSAGE: Directory', path, 'created.' |
---|
| 112 | |
---|
| 113 | except: |
---|
| 114 | print 'WARNING: Directory', path, 'could not be created.' |
---|
| 115 | if unix: |
---|
| 116 | path = '/tmp/' |
---|
| 117 | else: |
---|
| 118 | path = 'C:' |
---|
| 119 | |
---|
| 120 | print 'Using directory %s instead' %path |
---|
| 121 | |
---|
| 122 | return(path) |
---|
| 123 | |
---|
| 124 | |
---|
| 125 | |
---|
| 126 | def del_dir(path): |
---|
| 127 | """Recursively delete directory path and all its contents |
---|
| 128 | """ |
---|
| 129 | |
---|
| 130 | import os |
---|
| 131 | |
---|
| 132 | if os.path.isdir(path): |
---|
| 133 | for file in os.listdir(path): |
---|
| 134 | X = os.path.join(path, file) |
---|
| 135 | |
---|
| 136 | |
---|
| 137 | if os.path.isdir(X) and not os.path.islink(X): |
---|
| 138 | del_dir(X) |
---|
| 139 | else: |
---|
| 140 | try: |
---|
| 141 | os.remove(X) |
---|
| 142 | except: |
---|
| 143 | print "Could not remove file %s" %X |
---|
| 144 | |
---|
| 145 | os.rmdir(path) |
---|
| 146 | |
---|
| 147 | |
---|
| 148 | |
---|
| 149 | def create_filename(datadir, filename, format, size=None, time=None): |
---|
| 150 | |
---|
| 151 | import os |
---|
| 152 | #from config import data_dir |
---|
| 153 | |
---|
| 154 | FN = check_dir(datadir) + filename |
---|
| 155 | |
---|
| 156 | if size is not None: |
---|
| 157 | FN += '_size%d' %size |
---|
| 158 | |
---|
| 159 | if time is not None: |
---|
| 160 | FN += '_time%.2f' %time |
---|
| 161 | |
---|
| 162 | FN += '.' + format |
---|
| 163 | return FN |
---|
| 164 | |
---|
| 165 | |
---|
| 166 | def get_files(datadir, filename, format, size): |
---|
| 167 | """Get all file (names) with gven name, size and format |
---|
| 168 | """ |
---|
| 169 | |
---|
| 170 | import glob |
---|
| 171 | |
---|
| 172 | import os |
---|
| 173 | #from config import data_dir |
---|
| 174 | |
---|
| 175 | dir = check_dir(datadir) |
---|
| 176 | |
---|
| 177 | pattern = dir + os.sep + filename + '_size=%d*.%s' %(size, format) |
---|
| 178 | return glob.glob(pattern) |
---|
| 179 | |
---|
| 180 | |
---|
| 181 | |
---|
| 182 | #Generic class for storing output to e.g. visualisation or checkpointing |
---|
| 183 | class Data_format: |
---|
| 184 | """Generic interface to data formats |
---|
| 185 | """ |
---|
| 186 | |
---|
| 187 | |
---|
| 188 | def __init__(self, domain, extension, mode = 'w'): |
---|
| 189 | assert mode in ['r', 'w', 'a'], '''Mode %s must be either:''' %mode +\ |
---|
| 190 | ''' 'w' (write)'''+\ |
---|
| 191 | ''' 'r' (read)''' +\ |
---|
| 192 | ''' 'a' (append)''' |
---|
| 193 | |
---|
| 194 | #Create filename |
---|
| 195 | #self.filename = create_filename(domain.get_datadir(), |
---|
| 196 | #domain.get_name(), extension, len(domain)) |
---|
| 197 | |
---|
| 198 | |
---|
| 199 | self.filename = create_filename(domain.get_datadir(), |
---|
| 200 | domain.get_name(), extension) |
---|
| 201 | |
---|
| 202 | #print 'F', self.filename |
---|
| 203 | self.timestep = 0 |
---|
| 204 | #I CHANGED THIS |
---|
| 205 | #self.number_of_volumes = len(domain) |
---|
| 206 | self.number_of_volumes = domain.number_of_elements |
---|
| 207 | self.domain = domain |
---|
| 208 | |
---|
| 209 | |
---|
| 210 | #FIXME: Should we have a general set_precision function? |
---|
| 211 | |
---|
| 212 | |
---|
| 213 | |
---|
| 214 | #Class for storing output to e.g. visualisation |
---|
| 215 | class Data_format_sww(Data_format): |
---|
| 216 | """Interface to native NetCDF format (.sww) for storing model output |
---|
| 217 | |
---|
| 218 | There are two kinds of data |
---|
| 219 | |
---|
| 220 | 1: Constant data: Vertex coordinates and field values. Stored once |
---|
| 221 | 2: Variable data: Conserved quantities. Stored once per timestep. |
---|
| 222 | |
---|
| 223 | All data is assumed to reside at vertex locations. |
---|
| 224 | """ |
---|
| 225 | |
---|
| 226 | |
---|
| 227 | def __init__(self, domain, mode = 'w',\ |
---|
| 228 | max_size = 2000000000, |
---|
| 229 | recursion = False): |
---|
| 230 | from Scientific.IO.NetCDF import NetCDFFile |
---|
| 231 | from Numeric import Int, Float, Float32 |
---|
| 232 | |
---|
| 233 | self.precision = Float32 #Use single precision |
---|
| 234 | if hasattr(domain, 'max_size'): |
---|
| 235 | self.max_size = domain.max_size #file size max is 2Gig |
---|
| 236 | else: |
---|
| 237 | self.max_size = max_size |
---|
| 238 | self.recursion = recursion |
---|
| 239 | self.mode = mode |
---|
| 240 | |
---|
| 241 | Data_format.__init__(self, domain, 'sww', mode) |
---|
| 242 | |
---|
| 243 | |
---|
| 244 | # NetCDF file definition |
---|
| 245 | fid = NetCDFFile(self.filename, mode) |
---|
| 246 | |
---|
| 247 | if mode == 'w': |
---|
| 248 | |
---|
| 249 | #Create new file |
---|
| 250 | fid.institution = 'Geoscience Australia' |
---|
| 251 | fid.description = 'Output from pyvolution suitable for plotting' |
---|
| 252 | |
---|
| 253 | if domain.smooth: |
---|
| 254 | fid.smoothing = 'Yes' |
---|
| 255 | else: |
---|
| 256 | fid.smoothing = 'No' |
---|
| 257 | |
---|
| 258 | fid.order = domain.default_order |
---|
| 259 | |
---|
| 260 | if hasattr(domain, 'texture'): |
---|
| 261 | fid.texture = domain.texture |
---|
| 262 | #else: |
---|
| 263 | # fid.texture = 'None' |
---|
| 264 | |
---|
| 265 | #Reference point |
---|
| 266 | #Start time in seconds since the epoch (midnight 1/1/1970) |
---|
| 267 | #FIXME: Use Georef |
---|
| 268 | fid.starttime = domain.starttime |
---|
| 269 | |
---|
| 270 | # dimension definitions |
---|
| 271 | fid.createDimension('number_of_volumes', self.number_of_volumes) |
---|
| 272 | #fid.createDimension('number_of_vertices', 3) |
---|
| 273 | fid.createDimension('number_of_vertices', 2) |
---|
| 274 | |
---|
| 275 | if domain.smooth is True: |
---|
| 276 | #i changed this |
---|
| 277 | #fid.createDimension('number_of_points', len(domain.vertexlist)) |
---|
| 278 | fid.createDimension('number_of_points', domain.number_of_elements+1) |
---|
| 279 | else: |
---|
| 280 | #fid.createDimension('number_of_points', 3*self.number_of_volumes) |
---|
| 281 | fid.createDimension('number_of_points', 2*self.number_of_volumes) |
---|
| 282 | |
---|
| 283 | fid.createDimension('number_of_timesteps', None) #extensible |
---|
| 284 | |
---|
| 285 | # variable definitions |
---|
| 286 | fid.createVariable('x', self.precision, ('number_of_points',)) |
---|
| 287 | #fid.createVariable('y', self.precision, ('number_of_points',)) |
---|
| 288 | fid.createVariable('elevation', self.precision, ('number_of_points',)) |
---|
| 289 | if domain.geo_reference is not None: |
---|
| 290 | domain.geo_reference.write_NetCDF(fid) |
---|
| 291 | |
---|
| 292 | #FIXME: Backwards compatibility |
---|
| 293 | fid.createVariable('z', self.precision, ('number_of_points',)) |
---|
| 294 | ################################# |
---|
| 295 | |
---|
| 296 | fid.createVariable('volumes', Int, ('number_of_volumes', |
---|
| 297 | 'number_of_vertices')) |
---|
| 298 | |
---|
| 299 | fid.createVariable('time', self.precision, |
---|
| 300 | ('number_of_timesteps',)) |
---|
| 301 | |
---|
| 302 | fid.createVariable('stage', self.precision, |
---|
| 303 | ('number_of_timesteps', |
---|
| 304 | 'number_of_points')) |
---|
| 305 | |
---|
| 306 | fid.createVariable('xmomentum', self.precision, |
---|
| 307 | ('number_of_timesteps', |
---|
| 308 | 'number_of_points')) |
---|
| 309 | |
---|
| 310 | #fid.createVariable('ymomentum', self.precision, |
---|
| 311 | # ('number_of_timesteps', |
---|
| 312 | # 'number_of_points')) |
---|
| 313 | |
---|
| 314 | #Close |
---|
| 315 | fid.close() |
---|
| 316 | |
---|
| 317 | |
---|
| 318 | def store_connectivity(self): |
---|
| 319 | """Specialisation of store_connectivity for net CDF format |
---|
| 320 | |
---|
| 321 | Writes x,y,z coordinates of triangles constituting |
---|
| 322 | the bed elevation. |
---|
| 323 | """ |
---|
| 324 | |
---|
| 325 | from Scientific.IO.NetCDF import NetCDFFile |
---|
| 326 | |
---|
| 327 | from Numeric import concatenate, Int |
---|
| 328 | |
---|
| 329 | domain = self.domain |
---|
| 330 | |
---|
| 331 | #Get NetCDF |
---|
| 332 | fid = NetCDFFile(self.filename, 'a') #Open existing file for append |
---|
| 333 | |
---|
| 334 | # Get the variables |
---|
| 335 | x = fid.variables['x'] |
---|
| 336 | #y = fid.variables['y'] |
---|
| 337 | z = fid.variables['elevation'] |
---|
| 338 | |
---|
| 339 | volumes = fid.variables['volumes'] |
---|
| 340 | |
---|
| 341 | # Get X, Y and bed elevation Z |
---|
| 342 | Q = domain.quantities['elevation'] |
---|
| 343 | #X,Y,Z,V = Q.get_vertex_values(xy=True, |
---|
| 344 | # precision = self.precision) |
---|
| 345 | X,Z,V = Q.get_vertex_values(x=True, |
---|
| 346 | precision = self.precision) |
---|
| 347 | |
---|
| 348 | x[:] = X.astype(self.precision) |
---|
| 349 | #y[:] = Y.astype(self.precision) |
---|
| 350 | z[:] = Z.astype(self.precision) |
---|
| 351 | |
---|
| 352 | #FIXME: Backwards compatibility |
---|
| 353 | z = fid.variables['z'] |
---|
| 354 | z[:] = Z.astype(self.precision) |
---|
| 355 | ################################ |
---|
| 356 | |
---|
| 357 | volumes[:] = V.astype(volumes.typecode()) |
---|
| 358 | |
---|
| 359 | #Close |
---|
| 360 | fid.close() |
---|
| 361 | |
---|
| 362 | def store_timestep(self, names): |
---|
| 363 | """Store time and named quantities to file |
---|
| 364 | """ |
---|
| 365 | from Scientific.IO.NetCDF import NetCDFFile |
---|
| 366 | import types |
---|
| 367 | from time import sleep |
---|
| 368 | from os import stat |
---|
| 369 | |
---|
| 370 | #minimum_allowed_depth = 0.001 |
---|
| 371 | minimum_allowed_depth = 0.0 #FIXME pass in or read from domain |
---|
| 372 | from Numeric import choose |
---|
| 373 | |
---|
| 374 | #Get NetCDF |
---|
| 375 | retries = 0 |
---|
| 376 | file_open = False |
---|
| 377 | while not file_open and retries < 10: |
---|
| 378 | try: |
---|
| 379 | fid = NetCDFFile(self.filename, 'a') #Open existing file |
---|
| 380 | except IOError: |
---|
| 381 | #This could happen if someone was reading the file. |
---|
| 382 | #In that case, wait a while and try again |
---|
| 383 | msg = 'Warning (store_timestep): File %s could not be opened'\ |
---|
| 384 | %self.filename |
---|
| 385 | msg += ' - trying step %s again' %self.domain.time |
---|
| 386 | print msg |
---|
| 387 | retries += 1 |
---|
| 388 | sleep(1) |
---|
| 389 | else: |
---|
| 390 | file_open = True |
---|
| 391 | |
---|
| 392 | if not file_open: |
---|
| 393 | msg = 'File %s could not be opened for append' %self.filename |
---|
| 394 | raise msg |
---|
| 395 | |
---|
| 396 | |
---|
| 397 | |
---|
| 398 | #Check to see if the file is already too big: |
---|
| 399 | time = fid.variables['time'] |
---|
| 400 | i = len(time)+1 |
---|
| 401 | file_size = stat(self.filename)[6] |
---|
| 402 | file_size_increase = file_size/i |
---|
| 403 | if file_size + file_size_increase > self.max_size*(2**self.recursion): |
---|
| 404 | #in order to get the file name and start time correct, |
---|
| 405 | #I change the domian.filename and domain.starttime. |
---|
| 406 | #This is the only way to do this without changing |
---|
| 407 | #other modules (I think). |
---|
| 408 | |
---|
| 409 | #write a filename addon that won't break swollens reader |
---|
| 410 | #(10.sww is bad) |
---|
| 411 | filename_ext = '_time_%s'%self.domain.time |
---|
| 412 | filename_ext = filename_ext.replace('.', '_') |
---|
| 413 | #remember the old filename, then give domain a |
---|
| 414 | #name with the extension |
---|
| 415 | old_domain_filename = self.domain.filename |
---|
| 416 | if not self.recursion: |
---|
| 417 | self.domain.filename = self.domain.filename+filename_ext |
---|
| 418 | |
---|
| 419 | #change the domain starttime to the current time |
---|
| 420 | old_domain_starttime = self.domain.starttime |
---|
| 421 | self.domain.starttime = self.domain.time |
---|
| 422 | |
---|
| 423 | #build a new data_structure. |
---|
| 424 | next_data_structure=\ |
---|
| 425 | Data_format_sww(self.domain, mode=self.mode,\ |
---|
| 426 | max_size = self.max_size,\ |
---|
| 427 | recursion = self.recursion+1) |
---|
| 428 | if not self.recursion: |
---|
| 429 | print ' file_size = %s'%file_size |
---|
| 430 | print ' saving file to %s'%next_data_structure.filename |
---|
| 431 | #set up the new data_structure |
---|
| 432 | self.domain.writer = next_data_structure |
---|
| 433 | |
---|
| 434 | #FIXME - could be cleaner to use domain.store_timestep etc. |
---|
| 435 | next_data_structure.store_connectivity() |
---|
| 436 | next_data_structure.store_timestep(names) |
---|
| 437 | fid.sync() |
---|
| 438 | fid.close() |
---|
| 439 | |
---|
| 440 | #restore the old starttime and filename |
---|
| 441 | self.domain.starttime = old_domain_starttime |
---|
| 442 | self.domain.filename = old_domain_filename |
---|
| 443 | else: |
---|
| 444 | self.recursion = False |
---|
| 445 | domain = self.domain |
---|
| 446 | |
---|
| 447 | # Get the variables |
---|
| 448 | time = fid.variables['time'] |
---|
| 449 | stage = fid.variables['stage'] |
---|
| 450 | xmomentum = fid.variables['xmomentum'] |
---|
| 451 | #ymomentum = fid.variables['ymomentum'] |
---|
| 452 | i = len(time) |
---|
| 453 | |
---|
| 454 | #Store time |
---|
| 455 | time[i] = self.domain.time |
---|
| 456 | |
---|
| 457 | |
---|
| 458 | if type(names) not in [types.ListType, types.TupleType]: |
---|
| 459 | names = [names] |
---|
| 460 | |
---|
| 461 | for name in names: |
---|
| 462 | # Get quantity |
---|
| 463 | Q = domain.quantities[name] |
---|
| 464 | #A,V = Q.get_vertex_values(xy = False, |
---|
| 465 | # precision = self.precision) |
---|
| 466 | A,V = Q.get_vertex_values(x = False, |
---|
| 467 | precision = self.precision) |
---|
| 468 | |
---|
| 469 | #FIXME: Make this general (see below) |
---|
| 470 | if name == 'stage': |
---|
| 471 | z = fid.variables['elevation'] |
---|
| 472 | #print z[:] |
---|
| 473 | #print A-z[:] |
---|
| 474 | A = choose( A-z[:] >= minimum_allowed_depth, (z[:], A)) |
---|
| 475 | stage[i,:] = A.astype(self.precision) |
---|
| 476 | elif name == 'xmomentum': |
---|
| 477 | xmomentum[i,:] = A.astype(self.precision) |
---|
| 478 | elif name == 'ymomentum': |
---|
| 479 | ymomentum[i,:] = A.astype(self.precision) |
---|
| 480 | |
---|
| 481 | #As in.... |
---|
| 482 | #eval( name + '[i,:] = A.astype(self.precision)' ) |
---|
| 483 | #FIXME: But we need a UNIT test for that before refactoring |
---|
| 484 | |
---|
| 485 | |
---|
| 486 | |
---|
| 487 | #Flush and close |
---|
| 488 | fid.sync() |
---|
| 489 | fid.close() |
---|
| 490 | |
---|
| 491 | |
---|
| 492 | |
---|
| 493 | #Class for handling checkpoints data |
---|
| 494 | class Data_format_cpt(Data_format): |
---|
| 495 | """Interface to native NetCDF format (.cpt) |
---|
| 496 | """ |
---|
| 497 | |
---|
| 498 | |
---|
| 499 | def __init__(self, domain, mode = 'w'): |
---|
| 500 | from Scientific.IO.NetCDF import NetCDFFile |
---|
| 501 | from Numeric import Int, Float, Float |
---|
| 502 | |
---|
| 503 | self.precision = Float #Use full precision |
---|
| 504 | |
---|
| 505 | Data_format.__init__(self, domain, 'sww', mode) |
---|
| 506 | |
---|
| 507 | |
---|
| 508 | # NetCDF file definition |
---|
| 509 | fid = NetCDFFile(self.filename, mode) |
---|
| 510 | |
---|
| 511 | if mode == 'w': |
---|
| 512 | #Create new file |
---|
| 513 | fid.institution = 'Geoscience Australia' |
---|
| 514 | fid.description = 'Checkpoint data' |
---|
| 515 | #fid.smooth = domain.smooth |
---|
| 516 | fid.order = domain.default_order |
---|
| 517 | |
---|
| 518 | # dimension definitions |
---|
| 519 | fid.createDimension('number_of_volumes', self.number_of_volumes) |
---|
| 520 | fid.createDimension('number_of_vertices', 3) |
---|
| 521 | |
---|
| 522 | #Store info at all vertices (no smoothing) |
---|
| 523 | fid.createDimension('number_of_points', 3*self.number_of_volumes) |
---|
| 524 | fid.createDimension('number_of_timesteps', None) #extensible |
---|
| 525 | |
---|
| 526 | # variable definitions |
---|
| 527 | |
---|
| 528 | #Mesh |
---|
| 529 | fid.createVariable('x', self.precision, ('number_of_points',)) |
---|
| 530 | fid.createVariable('y', self.precision, ('number_of_points',)) |
---|
| 531 | |
---|
| 532 | |
---|
| 533 | fid.createVariable('volumes', Int, ('number_of_volumes', |
---|
| 534 | 'number_of_vertices')) |
---|
| 535 | |
---|
| 536 | fid.createVariable('time', self.precision, |
---|
| 537 | ('number_of_timesteps',)) |
---|
| 538 | |
---|
| 539 | #Allocate space for all quantities |
---|
| 540 | for name in domain.quantities.keys(): |
---|
| 541 | fid.createVariable(name, self.precision, |
---|
| 542 | ('number_of_timesteps', |
---|
| 543 | 'number_of_points')) |
---|
| 544 | |
---|
| 545 | #Close |
---|
| 546 | fid.close() |
---|
| 547 | |
---|
| 548 | |
---|
| 549 | def store_checkpoint(self): |
---|
| 550 | """ |
---|
| 551 | Write x,y coordinates of triangles. |
---|
| 552 | Write connectivity ( |
---|
| 553 | constituting |
---|
| 554 | the bed elevation. |
---|
| 555 | """ |
---|
| 556 | |
---|
| 557 | from Scientific.IO.NetCDF import NetCDFFile |
---|
| 558 | |
---|
| 559 | from Numeric import concatenate |
---|
| 560 | |
---|
| 561 | domain = self.domain |
---|
| 562 | |
---|
| 563 | #Get NetCDF |
---|
| 564 | fid = NetCDFFile(self.filename, 'a') #Open existing file for append |
---|
| 565 | |
---|
| 566 | # Get the variables |
---|
| 567 | x = fid.variables['x'] |
---|
| 568 | y = fid.variables['y'] |
---|
| 569 | |
---|
| 570 | volumes = fid.variables['volumes'] |
---|
| 571 | |
---|
| 572 | # Get X, Y and bed elevation Z |
---|
| 573 | Q = domain.quantities['elevation'] |
---|
| 574 | X,Y,Z,V = Q.get_vertex_values(xy=True, |
---|
| 575 | precision = self.precision) |
---|
| 576 | |
---|
| 577 | |
---|
| 578 | |
---|
| 579 | x[:] = X.astype(self.precision) |
---|
| 580 | y[:] = Y.astype(self.precision) |
---|
| 581 | z[:] = Z.astype(self.precision) |
---|
| 582 | |
---|
| 583 | volumes[:] = V |
---|
| 584 | |
---|
| 585 | #Close |
---|
| 586 | fid.close() |
---|
| 587 | |
---|
| 588 | |
---|
| 589 | def store_timestep(self, name): |
---|
| 590 | """Store time and named quantity to file |
---|
| 591 | """ |
---|
| 592 | from Scientific.IO.NetCDF import NetCDFFile |
---|
| 593 | from time import sleep |
---|
| 594 | |
---|
| 595 | #Get NetCDF |
---|
| 596 | retries = 0 |
---|
| 597 | file_open = False |
---|
| 598 | while not file_open and retries < 10: |
---|
| 599 | try: |
---|
| 600 | fid = NetCDFFile(self.filename, 'a') #Open existing file |
---|
| 601 | except IOError: |
---|
| 602 | #This could happen if someone was reading the file. |
---|
| 603 | #In that case, wait a while and try again |
---|
| 604 | msg = 'Warning (store_timestep): File %s could not be opened'\ |
---|
| 605 | %self.filename |
---|
| 606 | msg += ' - trying again' |
---|
| 607 | print msg |
---|
| 608 | retries += 1 |
---|
| 609 | sleep(1) |
---|
| 610 | else: |
---|
| 611 | file_open = True |
---|
| 612 | |
---|
| 613 | if not file_open: |
---|
| 614 | msg = 'File %s could not be opened for append' %self.filename |
---|
| 615 | raise msg |
---|
| 616 | |
---|
| 617 | |
---|
| 618 | domain = self.domain |
---|
| 619 | |
---|
| 620 | # Get the variables |
---|
| 621 | time = fid.variables['time'] |
---|
| 622 | stage = fid.variables['stage'] |
---|
| 623 | i = len(time) |
---|
| 624 | |
---|
| 625 | #Store stage |
---|
| 626 | time[i] = self.domain.time |
---|
| 627 | |
---|
| 628 | # Get quantity |
---|
| 629 | Q = domain.quantities[name] |
---|
| 630 | A,V = Q.get_vertex_values(xy=False, |
---|
| 631 | precision = self.precision) |
---|
| 632 | |
---|
| 633 | stage[i,:] = A.astype(self.precision) |
---|
| 634 | |
---|
| 635 | #Flush and close |
---|
| 636 | fid.sync() |
---|
| 637 | fid.close() |
---|
| 638 | |
---|
| 639 | |
---|
| 640 | |
---|
| 641 | |
---|
| 642 | |
---|
| 643 | #Function for storing xya output |
---|
| 644 | #FIXME Not done yet for this version |
---|
| 645 | class Data_format_xya(Data_format): |
---|
| 646 | """Generic interface to data formats |
---|
| 647 | """ |
---|
| 648 | |
---|
| 649 | |
---|
| 650 | def __init__(self, domain, mode = 'w'): |
---|
| 651 | from Scientific.IO.NetCDF import NetCDFFile |
---|
| 652 | from Numeric import Int, Float, Float32 |
---|
| 653 | |
---|
| 654 | self.precision = Float32 #Use single precision |
---|
| 655 | |
---|
| 656 | Data_format.__init__(self, domain, 'xya', mode) |
---|
| 657 | |
---|
| 658 | |
---|
| 659 | |
---|
| 660 | #FIXME -This is the old xya format |
---|
| 661 | def store_all(self): |
---|
| 662 | """Specialisation of store all for xya format |
---|
| 663 | |
---|
| 664 | Writes x,y,z coordinates of triangles constituting |
---|
| 665 | the bed elevation. |
---|
| 666 | """ |
---|
| 667 | |
---|
| 668 | from Numeric import concatenate |
---|
| 669 | |
---|
| 670 | domain = self.domain |
---|
| 671 | |
---|
| 672 | fd = open(self.filename, 'w') |
---|
| 673 | |
---|
| 674 | |
---|
| 675 | if domain.smooth is True: |
---|
| 676 | number_of_points = len(domain.vertexlist) |
---|
| 677 | else: |
---|
| 678 | number_of_points = 3*self.number_of_volumes |
---|
| 679 | |
---|
| 680 | numVertAttrib = 3 #Three attributes is what is assumed by the xya format |
---|
| 681 | |
---|
| 682 | fd.write(str(number_of_points) + " " + str(numVertAttrib) +\ |
---|
| 683 | " # <vertex #> <x> <y> [attributes]" + "\n") |
---|
| 684 | |
---|
| 685 | |
---|
| 686 | # Get X, Y, bed elevation and friction (index=0,1) |
---|
| 687 | X,Y,A,V = domain.get_vertex_values(xy=True, value_array='field_values', |
---|
| 688 | indices = (0,1), precision = self.precision) |
---|
| 689 | |
---|
| 690 | bed_eles = A[:,0] |
---|
| 691 | fricts = A[:,1] |
---|
| 692 | |
---|
| 693 | # Get stage (index=0) |
---|
| 694 | B,V = domain.get_vertex_values(xy=False, value_array='conserved_quantities', |
---|
| 695 | indices = (0,), precision = self.precision) |
---|
| 696 | |
---|
| 697 | stages = B[:,0] |
---|
| 698 | |
---|
| 699 | #<vertex #> <x> <y> [attributes] |
---|
| 700 | for x, y, bed_ele, stage, frict in map(None, X, Y, bed_eles, |
---|
| 701 | stages, fricts): |
---|
| 702 | |
---|
| 703 | s = '%.6f %.6f %.6f %.6f %.6f\n' %(x, y, bed_ele, stage, frict) |
---|
| 704 | fd.write(s) |
---|
| 705 | |
---|
| 706 | #close |
---|
| 707 | fd.close() |
---|
| 708 | |
---|
| 709 | |
---|
| 710 | def store_timestep(self, t, V0, V1, V2): |
---|
| 711 | """Store time, water heights (and momentums) to file |
---|
| 712 | """ |
---|
| 713 | pass |
---|
| 714 | |
---|
| 715 | |
---|
| 716 | #Auxiliary |
---|
| 717 | def write_obj(filename,x,y,z): |
---|
| 718 | """Store x,y,z vectors into filename (obj format) |
---|
| 719 | Vectors are assumed to have dimension (M,3) where |
---|
| 720 | M corresponds to the number elements. |
---|
| 721 | triangles are assumed to be disconnected |
---|
| 722 | |
---|
| 723 | The three numbers in each vector correspond to three vertices, |
---|
| 724 | |
---|
| 725 | e.g. the x coordinate of vertex 1 of element i is in x[i,1] |
---|
| 726 | |
---|
| 727 | """ |
---|
| 728 | #print 'Writing obj to %s' % filename |
---|
| 729 | |
---|
| 730 | import os.path |
---|
| 731 | |
---|
| 732 | root, ext = os.path.splitext(filename) |
---|
| 733 | if ext == '.obj': |
---|
| 734 | FN = filename |
---|
| 735 | else: |
---|
| 736 | FN = filename + '.obj' |
---|
| 737 | |
---|
| 738 | |
---|
| 739 | outfile = open(FN, 'wb') |
---|
| 740 | outfile.write("# Triangulation as an obj file\n") |
---|
| 741 | |
---|
| 742 | M, N = x.shape |
---|
| 743 | assert N==3 #Assuming three vertices per element |
---|
| 744 | |
---|
| 745 | for i in range(M): |
---|
| 746 | for j in range(N): |
---|
| 747 | outfile.write("v %f %f %f\n" % (x[i,j],y[i,j],z[i,j])) |
---|
| 748 | |
---|
| 749 | for i in range(M): |
---|
| 750 | base = i*N |
---|
| 751 | outfile.write("f %d %d %d\n" % (base+1,base+2,base+3)) |
---|
| 752 | |
---|
| 753 | outfile.close() |
---|
| 754 | |
---|
| 755 | |
---|
| 756 | ######################################################### |
---|
| 757 | #Conversion routines |
---|
| 758 | ######################################################## |
---|
| 759 | |
---|
| 760 | def sww2obj(basefilename, size): |
---|
| 761 | """Convert netcdf based data output to obj |
---|
| 762 | """ |
---|
| 763 | from Scientific.IO.NetCDF import NetCDFFile |
---|
| 764 | |
---|
| 765 | from Numeric import Float, zeros |
---|
| 766 | |
---|
| 767 | #Get NetCDF |
---|
| 768 | FN = create_filename('.', basefilename, 'sww', size) |
---|
| 769 | print 'Reading from ', FN |
---|
| 770 | fid = NetCDFFile(FN, 'r') #Open existing file for read |
---|
| 771 | |
---|
| 772 | |
---|
| 773 | # Get the variables |
---|
| 774 | x = fid.variables['x'] |
---|
| 775 | y = fid.variables['y'] |
---|
| 776 | z = fid.variables['elevation'] |
---|
| 777 | time = fid.variables['time'] |
---|
| 778 | stage = fid.variables['stage'] |
---|
| 779 | |
---|
| 780 | M = size #Number of lines |
---|
| 781 | xx = zeros((M,3), Float) |
---|
| 782 | yy = zeros((M,3), Float) |
---|
| 783 | zz = zeros((M,3), Float) |
---|
| 784 | |
---|
| 785 | for i in range(M): |
---|
| 786 | for j in range(3): |
---|
| 787 | xx[i,j] = x[i+j*M] |
---|
| 788 | yy[i,j] = y[i+j*M] |
---|
| 789 | zz[i,j] = z[i+j*M] |
---|
| 790 | |
---|
| 791 | #Write obj for bathymetry |
---|
| 792 | FN = create_filename('.', basefilename, 'obj', size) |
---|
| 793 | write_obj(FN,xx,yy,zz) |
---|
| 794 | |
---|
| 795 | |
---|
| 796 | #Now read all the data with variable information, combine with |
---|
| 797 | #x,y info and store as obj |
---|
| 798 | |
---|
| 799 | for k in range(len(time)): |
---|
| 800 | t = time[k] |
---|
| 801 | print 'Processing timestep %f' %t |
---|
| 802 | |
---|
| 803 | for i in range(M): |
---|
| 804 | for j in range(3): |
---|
| 805 | zz[i,j] = stage[k,i+j*M] |
---|
| 806 | |
---|
| 807 | |
---|
| 808 | #Write obj for variable data |
---|
| 809 | #FN = create_filename(basefilename, 'obj', size, time=t) |
---|
| 810 | FN = create_filename('.', basefilename[:5], 'obj', size, time=t) |
---|
| 811 | write_obj(FN,xx,yy,zz) |
---|
| 812 | |
---|
| 813 | |
---|
| 814 | def dat2obj(basefilename): |
---|
| 815 | """Convert line based data output to obj |
---|
| 816 | FIXME: Obsolete? |
---|
| 817 | """ |
---|
| 818 | |
---|
| 819 | import glob, os |
---|
| 820 | from config import data_dir |
---|
| 821 | |
---|
| 822 | |
---|
| 823 | #Get bathymetry and x,y's |
---|
| 824 | lines = open(data_dir+os.sep+basefilename+'_geometry.dat', 'r').readlines() |
---|
| 825 | |
---|
| 826 | from Numeric import zeros, Float |
---|
| 827 | |
---|
| 828 | M = len(lines) #Number of lines |
---|
| 829 | x = zeros((M,3), Float) |
---|
| 830 | y = zeros((M,3), Float) |
---|
| 831 | z = zeros((M,3), Float) |
---|
| 832 | |
---|
| 833 | ##i = 0 |
---|
| 834 | for i, line in enumerate(lines): |
---|
| 835 | tokens = line.split() |
---|
| 836 | values = map(float,tokens) |
---|
| 837 | |
---|
| 838 | for j in range(3): |
---|
| 839 | x[i,j] = values[j*3] |
---|
| 840 | y[i,j] = values[j*3+1] |
---|
| 841 | z[i,j] = values[j*3+2] |
---|
| 842 | |
---|
| 843 | ##i += 1 |
---|
| 844 | |
---|
| 845 | |
---|
| 846 | #Write obj for bathymetry |
---|
| 847 | write_obj(data_dir+os.sep+basefilename+'_geometry',x,y,z) |
---|
| 848 | |
---|
| 849 | |
---|
| 850 | #Now read all the data files with variable information, combine with |
---|
| 851 | #x,y info |
---|
| 852 | #and store as obj |
---|
| 853 | |
---|
| 854 | files = glob.glob(data_dir+os.sep+basefilename+'*.dat') |
---|
| 855 | |
---|
| 856 | for filename in files: |
---|
| 857 | print 'Processing %s' % filename |
---|
| 858 | |
---|
| 859 | lines = open(data_dir+os.sep+filename,'r').readlines() |
---|
| 860 | assert len(lines) == M |
---|
| 861 | root, ext = os.path.splitext(filename) |
---|
| 862 | |
---|
| 863 | #Get time from filename |
---|
| 864 | i0 = filename.find('_time=') |
---|
| 865 | if i0 == -1: |
---|
| 866 | #Skip bathymetry file |
---|
| 867 | continue |
---|
| 868 | |
---|
| 869 | i0 += 6 #Position where time starts |
---|
| 870 | i1 = filename.find('.dat') |
---|
| 871 | |
---|
| 872 | if i1 > i0: |
---|
| 873 | t = float(filename[i0:i1]) |
---|
| 874 | else: |
---|
| 875 | raise 'Hmmmm' |
---|
| 876 | |
---|
| 877 | |
---|
| 878 | |
---|
| 879 | ##i = 0 |
---|
| 880 | for i, line in enumerate(lines): |
---|
| 881 | tokens = line.split() |
---|
| 882 | values = map(float,tokens) |
---|
| 883 | |
---|
| 884 | for j in range(3): |
---|
| 885 | z[i,j] = values[j] |
---|
| 886 | |
---|
| 887 | ##i += 1 |
---|
| 888 | |
---|
| 889 | #Write obj for variable data |
---|
| 890 | write_obj(data_dir+os.sep+basefilename+'_time=%.4f' %t,x,y,z) |
---|
| 891 | |
---|
| 892 | |
---|
| 893 | def filter_netcdf(filename1, filename2, first=0, last=None, step = 1): |
---|
| 894 | """Read netcdf filename1, pick timesteps first:step:last and save to |
---|
| 895 | nettcdf file filename2 |
---|
| 896 | """ |
---|
| 897 | from Scientific.IO.NetCDF import NetCDFFile |
---|
| 898 | |
---|
| 899 | #Get NetCDF |
---|
| 900 | infile = NetCDFFile(filename1, 'r') #Open existing file for read |
---|
| 901 | outfile = NetCDFFile(filename2, 'w') #Open new file |
---|
| 902 | |
---|
| 903 | |
---|
| 904 | #Copy dimensions |
---|
| 905 | for d in infile.dimensions: |
---|
| 906 | outfile.createDimension(d, infile.dimensions[d]) |
---|
| 907 | |
---|
| 908 | for name in infile.variables: |
---|
| 909 | var = infile.variables[name] |
---|
| 910 | outfile.createVariable(name, var.typecode(), var.dimensions) |
---|
| 911 | |
---|
| 912 | |
---|
| 913 | #Copy the static variables |
---|
| 914 | for name in infile.variables: |
---|
| 915 | if name == 'time' or name == 'stage': |
---|
| 916 | pass |
---|
| 917 | else: |
---|
| 918 | #Copy |
---|
| 919 | outfile.variables[name][:] = infile.variables[name][:] |
---|
| 920 | |
---|
| 921 | #Copy selected timesteps |
---|
| 922 | time = infile.variables['time'] |
---|
| 923 | stage = infile.variables['stage'] |
---|
| 924 | |
---|
| 925 | newtime = outfile.variables['time'] |
---|
| 926 | newstage = outfile.variables['stage'] |
---|
| 927 | |
---|
| 928 | if last is None: |
---|
| 929 | last = len(time) |
---|
| 930 | |
---|
| 931 | selection = range(first, last, step) |
---|
| 932 | for i, j in enumerate(selection): |
---|
| 933 | print 'Copying timestep %d of %d (%f)' %(j, last-first, time[j]) |
---|
| 934 | newtime[i] = time[j] |
---|
| 935 | newstage[i,:] = stage[j,:] |
---|
| 936 | |
---|
| 937 | #Close |
---|
| 938 | infile.close() |
---|
| 939 | outfile.close() |
---|
| 940 | |
---|
| 941 | |
---|
| 942 | #Get data objects |
---|
| 943 | def get_dataobject(domain, mode='w'): |
---|
| 944 | """Return instance of class of given format using filename |
---|
| 945 | """ |
---|
| 946 | |
---|
| 947 | cls = eval('Data_format_%s' %domain.format) |
---|
| 948 | return cls(domain, mode) |
---|
| 949 | |
---|
| 950 | def xya2pts(basename_in, basename_out=None, verbose=False, |
---|
| 951 | #easting_min=None, easting_max=None, |
---|
| 952 | #northing_min=None, northing_max=None, |
---|
| 953 | stride = 1, |
---|
| 954 | attribute_name = 'elevation', |
---|
| 955 | z_func = None): |
---|
| 956 | """Read points data from ascii (.xya) |
---|
| 957 | |
---|
| 958 | Example: |
---|
| 959 | |
---|
| 960 | x(m) y(m) z(m) |
---|
| 961 | 0.00000e+00 0.00000e+00 1.3535000e-01 |
---|
| 962 | 0.00000e+00 1.40000e-02 1.3535000e-01 |
---|
| 963 | |
---|
| 964 | |
---|
| 965 | |
---|
| 966 | Convert to NetCDF pts format which is |
---|
| 967 | |
---|
| 968 | points: (Nx2) Float array |
---|
| 969 | elevation: N Float array |
---|
| 970 | |
---|
| 971 | Only lines that contain three numeric values are processed |
---|
| 972 | |
---|
| 973 | If z_func is specified, it will be applied to the third field |
---|
| 974 | """ |
---|
| 975 | |
---|
| 976 | import os |
---|
| 977 | #from Scientific.IO.NetCDF import NetCDFFile |
---|
| 978 | from Numeric import Float, arrayrange, concatenate |
---|
| 979 | |
---|
| 980 | root, ext = os.path.splitext(basename_in) |
---|
| 981 | |
---|
| 982 | if ext == '': ext = '.xya' |
---|
| 983 | |
---|
| 984 | #Get NetCDF |
---|
| 985 | infile = open(root + ext, 'r') #Open existing xya file for read |
---|
| 986 | |
---|
| 987 | if verbose: print 'Reading xya points from %s' %(root + ext) |
---|
| 988 | |
---|
| 989 | points = [] |
---|
| 990 | attribute = [] |
---|
| 991 | for i, line in enumerate(infile.readlines()): |
---|
| 992 | |
---|
| 993 | if i % stride != 0: continue |
---|
| 994 | |
---|
| 995 | fields = line.split() |
---|
| 996 | |
---|
| 997 | try: |
---|
| 998 | assert len(fields) == 3 |
---|
| 999 | except: |
---|
| 1000 | print 'WARNING: Line %d doesn\'t have 3 elements: %s' %(i, line) |
---|
| 1001 | |
---|
| 1002 | try: |
---|
| 1003 | x = float( fields[0] ) |
---|
| 1004 | y = float( fields[1] ) |
---|
| 1005 | z = float( fields[2] ) |
---|
| 1006 | except: |
---|
| 1007 | continue |
---|
| 1008 | |
---|
| 1009 | points.append( [x, y] ) |
---|
| 1010 | |
---|
| 1011 | if callable(z_func): |
---|
| 1012 | attribute.append(z_func(z)) |
---|
| 1013 | else: |
---|
| 1014 | attribute.append(z) |
---|
| 1015 | |
---|
| 1016 | |
---|
| 1017 | #Get output file |
---|
| 1018 | if basename_out == None: |
---|
| 1019 | ptsname = root + '.pts' |
---|
| 1020 | else: |
---|
| 1021 | ptsname = basename_out + '.pts' |
---|
| 1022 | |
---|
| 1023 | if verbose: print 'Store to NetCDF file %s' %ptsname |
---|
| 1024 | write_ptsfile(ptsname, points, attribute, attribute_name) |
---|
| 1025 | |
---|
| 1026 | |
---|
| 1027 | |
---|
| 1028 | ######Obsoleted by export_points in load_mesh |
---|
| 1029 | def write_ptsfile(ptsname, points, attribute, attribute_name = None, |
---|
| 1030 | zone=None, xllcorner=None, yllcorner=None): |
---|
| 1031 | """Write points and associated attribute to pts (NetCDF) format |
---|
| 1032 | """ |
---|
| 1033 | |
---|
| 1034 | print 'WARNING: write_ptsfile is obsolete. Use export_points from load_mesh.loadASCII instead.' |
---|
| 1035 | |
---|
| 1036 | from Numeric import Float |
---|
| 1037 | |
---|
| 1038 | if attribute_name is None: |
---|
| 1039 | attribute_name = 'attribute' |
---|
| 1040 | |
---|
| 1041 | |
---|
| 1042 | from Scientific.IO.NetCDF import NetCDFFile |
---|
| 1043 | |
---|
| 1044 | # NetCDF file definition |
---|
| 1045 | outfile = NetCDFFile(ptsname, 'w') |
---|
| 1046 | |
---|
| 1047 | |
---|
| 1048 | #Create new file |
---|
| 1049 | outfile.institution = 'Geoscience Australia' |
---|
| 1050 | outfile.description = 'NetCDF pts format for compact and '\ |
---|
| 1051 | 'portable storage of spatial point data' |
---|
| 1052 | |
---|
| 1053 | |
---|
| 1054 | #Georeferencing |
---|
| 1055 | from coordinate_transforms.geo_reference import Geo_reference |
---|
| 1056 | if zone is None: |
---|
| 1057 | assert xllcorner is None, 'xllcorner must be None' |
---|
| 1058 | assert yllcorner is None, 'yllcorner must be None' |
---|
| 1059 | Geo_reference().write_NetCDF(outfile) |
---|
| 1060 | else: |
---|
| 1061 | Geo_reference(zone, xllcorner, yllcorner).write_NetCDF(outfile) |
---|
| 1062 | |
---|
| 1063 | |
---|
| 1064 | |
---|
| 1065 | outfile.createDimension('number_of_points', len(points)) |
---|
| 1066 | outfile.createDimension('number_of_dimensions', 2) #This is 2d data |
---|
| 1067 | |
---|
| 1068 | # variable definitions |
---|
| 1069 | outfile.createVariable('points', Float, ('number_of_points', |
---|
| 1070 | 'number_of_dimensions')) |
---|
| 1071 | outfile.createVariable(attribute_name, Float, ('number_of_points',)) |
---|
| 1072 | |
---|
| 1073 | # Get handles to the variables |
---|
| 1074 | nc_points = outfile.variables['points'] |
---|
| 1075 | nc_attribute = outfile.variables[attribute_name] |
---|
| 1076 | |
---|
| 1077 | #Store data |
---|
| 1078 | nc_points[:, :] = points |
---|
| 1079 | nc_attribute[:] = attribute |
---|
| 1080 | |
---|
| 1081 | outfile.close() |
---|
| 1082 | |
---|
| 1083 | |
---|
| 1084 | def dem2pts(basename_in, basename_out=None, verbose=False, |
---|
| 1085 | easting_min=None, easting_max=None, |
---|
| 1086 | northing_min=None, northing_max=None): |
---|
| 1087 | """Read Digitial Elevation model from the following NetCDF format (.dem) |
---|
| 1088 | |
---|
| 1089 | Example: |
---|
| 1090 | |
---|
| 1091 | ncols 3121 |
---|
| 1092 | nrows 1800 |
---|
| 1093 | xllcorner 722000 |
---|
| 1094 | yllcorner 5893000 |
---|
| 1095 | cellsize 25 |
---|
| 1096 | NODATA_value -9999 |
---|
| 1097 | 138.3698 137.4194 136.5062 135.5558 .......... |
---|
| 1098 | |
---|
| 1099 | Convert to NetCDF pts format which is |
---|
| 1100 | |
---|
| 1101 | points: (Nx2) Float array |
---|
| 1102 | elevation: N Float array |
---|
| 1103 | """ |
---|
| 1104 | |
---|
| 1105 | #FIXME: Can this be written feasibly using write_pts? |
---|
| 1106 | |
---|
| 1107 | import os |
---|
| 1108 | from Scientific.IO.NetCDF import NetCDFFile |
---|
| 1109 | from Numeric import Float, zeros, reshape |
---|
| 1110 | |
---|
| 1111 | root = basename_in |
---|
| 1112 | |
---|
| 1113 | #Get NetCDF |
---|
| 1114 | infile = NetCDFFile(root + '.dem', 'r') #Open existing netcdf file for read |
---|
| 1115 | |
---|
| 1116 | if verbose: print 'Reading DEM from %s' %(root + '.dem') |
---|
| 1117 | |
---|
| 1118 | ncols = infile.ncols[0] |
---|
| 1119 | nrows = infile.nrows[0] |
---|
| 1120 | xllcorner = infile.xllcorner[0] #Easting of lower left corner |
---|
| 1121 | yllcorner = infile.yllcorner[0] #Northing of lower left corner |
---|
| 1122 | cellsize = infile.cellsize[0] |
---|
| 1123 | NODATA_value = infile.NODATA_value[0] |
---|
| 1124 | dem_elevation = infile.variables['elevation'] |
---|
| 1125 | |
---|
| 1126 | zone = infile.zone[0] |
---|
| 1127 | false_easting = infile.false_easting[0] |
---|
| 1128 | false_northing = infile.false_northing[0] |
---|
| 1129 | |
---|
| 1130 | #Text strings |
---|
| 1131 | projection = infile.projection |
---|
| 1132 | datum = infile.datum |
---|
| 1133 | units = infile.units |
---|
| 1134 | |
---|
| 1135 | |
---|
| 1136 | #Get output file |
---|
| 1137 | if basename_out == None: |
---|
| 1138 | ptsname = root + '.pts' |
---|
| 1139 | else: |
---|
| 1140 | ptsname = basename_out + '.pts' |
---|
| 1141 | |
---|
| 1142 | if verbose: print 'Store to NetCDF file %s' %ptsname |
---|
| 1143 | # NetCDF file definition |
---|
| 1144 | outfile = NetCDFFile(ptsname, 'w') |
---|
| 1145 | |
---|
| 1146 | #Create new file |
---|
| 1147 | outfile.institution = 'Geoscience Australia' |
---|
| 1148 | outfile.description = 'NetCDF pts format for compact and portable storage ' +\ |
---|
| 1149 | 'of spatial point data' |
---|
| 1150 | #assign default values |
---|
| 1151 | if easting_min is None: easting_min = xllcorner |
---|
| 1152 | if easting_max is None: easting_max = xllcorner + ncols*cellsize |
---|
| 1153 | if northing_min is None: northing_min = yllcorner |
---|
| 1154 | if northing_max is None: northing_max = yllcorner + nrows*cellsize |
---|
| 1155 | |
---|
| 1156 | #compute offsets to update georeferencing |
---|
| 1157 | easting_offset = xllcorner - easting_min |
---|
| 1158 | northing_offset = yllcorner - northing_min |
---|
| 1159 | |
---|
| 1160 | #Georeferencing |
---|
| 1161 | outfile.zone = zone |
---|
| 1162 | outfile.xllcorner = easting_min #Easting of lower left corner |
---|
| 1163 | outfile.yllcorner = northing_min #Northing of lower left corner |
---|
| 1164 | outfile.false_easting = false_easting |
---|
| 1165 | outfile.false_northing = false_northing |
---|
| 1166 | |
---|
| 1167 | outfile.projection = projection |
---|
| 1168 | outfile.datum = datum |
---|
| 1169 | outfile.units = units |
---|
| 1170 | |
---|
| 1171 | |
---|
| 1172 | #Grid info (FIXME: probably not going to be used, but heck) |
---|
| 1173 | outfile.ncols = ncols |
---|
| 1174 | outfile.nrows = nrows |
---|
| 1175 | |
---|
| 1176 | |
---|
| 1177 | # dimension definitions |
---|
| 1178 | nrows_in_bounding_box = int(round((northing_max-northing_min)/cellsize)) |
---|
| 1179 | ncols_in_bounding_box = int(round((easting_max-easting_min)/cellsize)) |
---|
| 1180 | outfile.createDimension('number_of_points', nrows_in_bounding_box*ncols_in_bounding_box) |
---|
| 1181 | outfile.createDimension('number_of_dimensions', 2) #This is 2d data |
---|
| 1182 | |
---|
| 1183 | # variable definitions |
---|
| 1184 | outfile.createVariable('points', Float, ('number_of_points', |
---|
| 1185 | 'number_of_dimensions')) |
---|
| 1186 | outfile.createVariable('elevation', Float, ('number_of_points',)) |
---|
| 1187 | |
---|
| 1188 | # Get handles to the variables |
---|
| 1189 | points = outfile.variables['points'] |
---|
| 1190 | elevation = outfile.variables['elevation'] |
---|
| 1191 | |
---|
| 1192 | dem_elevation_r = reshape(dem_elevation, (nrows, ncols)) |
---|
| 1193 | |
---|
| 1194 | #Store data |
---|
| 1195 | #FIXME: Could perhaps be faster using array operations (Fixed 27/7/05) |
---|
| 1196 | global_index = 0 |
---|
| 1197 | for i in range(nrows): |
---|
| 1198 | if verbose and i%((nrows+10)/10)==0: |
---|
| 1199 | print 'Processing row %d of %d' %(i, nrows) |
---|
| 1200 | |
---|
| 1201 | lower_index = global_index |
---|
| 1202 | tpoints = zeros((ncols_in_bounding_box, 2), Float) |
---|
| 1203 | telev = zeros(ncols_in_bounding_box, Float) |
---|
| 1204 | local_index = 0 |
---|
| 1205 | |
---|
| 1206 | y = (nrows-i)*cellsize + yllcorner |
---|
| 1207 | for j in range(ncols): |
---|
| 1208 | |
---|
| 1209 | x = j*cellsize + xllcorner |
---|
| 1210 | if easting_min <= x <= easting_max and \ |
---|
| 1211 | northing_min <= y <= northing_max: |
---|
| 1212 | tpoints[local_index, :] = [x-easting_min,y-northing_min] |
---|
| 1213 | telev[local_index] = dem_elevation_r[i, j] |
---|
| 1214 | global_index += 1 |
---|
| 1215 | local_index += 1 |
---|
| 1216 | |
---|
| 1217 | upper_index = global_index |
---|
| 1218 | |
---|
| 1219 | if upper_index == lower_index + ncols_in_bounding_box: |
---|
| 1220 | points[lower_index:upper_index, :] = tpoints |
---|
| 1221 | elevation[lower_index:upper_index] = telev |
---|
| 1222 | |
---|
| 1223 | assert global_index == nrows_in_bounding_box*ncols_in_bounding_box, 'index not equal to number of points' |
---|
| 1224 | |
---|
| 1225 | infile.close() |
---|
| 1226 | outfile.close() |
---|
| 1227 | |
---|
| 1228 | |
---|
| 1229 | |
---|
| 1230 | def _read_hecras_cross_sections(lines): |
---|
| 1231 | """Return block of surface lines for each cross section |
---|
| 1232 | Starts with SURFACE LINE, |
---|
| 1233 | Ends with END CROSS-SECTION |
---|
| 1234 | """ |
---|
| 1235 | |
---|
| 1236 | points = [] |
---|
| 1237 | |
---|
| 1238 | reading_surface = False |
---|
| 1239 | for i, line in enumerate(lines): |
---|
| 1240 | |
---|
| 1241 | if len(line.strip()) == 0: #Ignore blanks |
---|
| 1242 | continue |
---|
| 1243 | |
---|
| 1244 | if lines[i].strip().startswith('SURFACE LINE'): |
---|
| 1245 | reading_surface = True |
---|
| 1246 | continue |
---|
| 1247 | |
---|
| 1248 | if lines[i].strip().startswith('END') and reading_surface: |
---|
| 1249 | yield points |
---|
| 1250 | reading_surface = False |
---|
| 1251 | points = [] |
---|
| 1252 | |
---|
| 1253 | if reading_surface: |
---|
| 1254 | fields = line.strip().split(',') |
---|
| 1255 | easting = float(fields[0]) |
---|
| 1256 | northing = float(fields[1]) |
---|
| 1257 | elevation = float(fields[2]) |
---|
| 1258 | points.append([easting, northing, elevation]) |
---|
| 1259 | |
---|
| 1260 | |
---|
| 1261 | |
---|
| 1262 | |
---|
| 1263 | def hecras_cross_sections2pts(basename_in, |
---|
| 1264 | basename_out=None, |
---|
| 1265 | verbose=False): |
---|
| 1266 | """Read HEC-RAS Elevation datal from the following ASCII format (.sdf) |
---|
| 1267 | |
---|
| 1268 | Example: |
---|
| 1269 | |
---|
| 1270 | |
---|
| 1271 | # RAS export file created on Mon 15Aug2005 11:42 |
---|
| 1272 | # by HEC-RAS Version 3.1.1 |
---|
| 1273 | |
---|
| 1274 | BEGIN HEADER: |
---|
| 1275 | UNITS: METRIC |
---|
| 1276 | DTM TYPE: TIN |
---|
| 1277 | DTM: v:\1\cit\perth_topo\river_tin |
---|
| 1278 | STREAM LAYER: c:\local\hecras\21_02_03\up_canning_cent3d.shp |
---|
| 1279 | CROSS-SECTION LAYER: c:\local\hecras\21_02_03\up_can_xs3d.shp |
---|
| 1280 | MAP PROJECTION: UTM |
---|
| 1281 | PROJECTION ZONE: 50 |
---|
| 1282 | DATUM: AGD66 |
---|
| 1283 | VERTICAL DATUM: |
---|
| 1284 | NUMBER OF REACHES: 19 |
---|
| 1285 | NUMBER OF CROSS-SECTIONS: 14206 |
---|
| 1286 | END HEADER: |
---|
| 1287 | |
---|
| 1288 | |
---|
| 1289 | Only the SURFACE LINE data of the following form will be utilised |
---|
| 1290 | |
---|
| 1291 | CROSS-SECTION: |
---|
| 1292 | STREAM ID:Southern-Wungong |
---|
| 1293 | REACH ID:Southern-Wungong |
---|
| 1294 | STATION:19040.* |
---|
| 1295 | CUT LINE: |
---|
| 1296 | 405548.671603161 , 6438142.7594925 |
---|
| 1297 | 405734.536092045 , 6438326.10404912 |
---|
| 1298 | 405745.130459356 , 6438331.48627354 |
---|
| 1299 | 405813.89633823 , 6438368.6272789 |
---|
| 1300 | SURFACE LINE: |
---|
| 1301 | 405548.67, 6438142.76, 35.37 |
---|
| 1302 | 405552.24, 6438146.28, 35.41 |
---|
| 1303 | 405554.78, 6438148.78, 35.44 |
---|
| 1304 | 405555.80, 6438149.79, 35.44 |
---|
| 1305 | 405559.37, 6438153.31, 35.45 |
---|
| 1306 | 405560.88, 6438154.81, 35.44 |
---|
| 1307 | 405562.93, 6438156.83, 35.42 |
---|
| 1308 | 405566.50, 6438160.35, 35.38 |
---|
| 1309 | 405566.99, 6438160.83, 35.37 |
---|
| 1310 | ... |
---|
| 1311 | END CROSS-SECTION |
---|
| 1312 | |
---|
| 1313 | Convert to NetCDF pts format which is |
---|
| 1314 | |
---|
| 1315 | points: (Nx2) Float array |
---|
| 1316 | elevation: N Float array |
---|
| 1317 | """ |
---|
| 1318 | |
---|
| 1319 | #FIXME: Can this be written feasibly using write_pts? |
---|
| 1320 | |
---|
| 1321 | import os |
---|
| 1322 | from Scientific.IO.NetCDF import NetCDFFile |
---|
| 1323 | from Numeric import Float, zeros, reshape |
---|
| 1324 | |
---|
| 1325 | root = basename_in |
---|
| 1326 | |
---|
| 1327 | #Get ASCII file |
---|
| 1328 | infile = open(root + '.sdf', 'r') #Open SDF file for read |
---|
| 1329 | |
---|
| 1330 | if verbose: print 'Reading DEM from %s' %(root + '.sdf') |
---|
| 1331 | |
---|
| 1332 | lines = infile.readlines() |
---|
| 1333 | infile.close() |
---|
| 1334 | |
---|
| 1335 | if verbose: print 'Converting to pts format' |
---|
| 1336 | |
---|
| 1337 | i = 0 |
---|
| 1338 | while lines[i].strip() == '' or lines[i].strip().startswith('#'): |
---|
| 1339 | i += 1 |
---|
| 1340 | |
---|
| 1341 | assert lines[i].strip().upper() == 'BEGIN HEADER:' |
---|
| 1342 | i += 1 |
---|
| 1343 | |
---|
| 1344 | assert lines[i].strip().upper().startswith('UNITS:') |
---|
| 1345 | units = lines[i].strip().split()[1] |
---|
| 1346 | i += 1 |
---|
| 1347 | |
---|
| 1348 | assert lines[i].strip().upper().startswith('DTM TYPE:') |
---|
| 1349 | i += 1 |
---|
| 1350 | |
---|
| 1351 | assert lines[i].strip().upper().startswith('DTM:') |
---|
| 1352 | i += 1 |
---|
| 1353 | |
---|
| 1354 | assert lines[i].strip().upper().startswith('STREAM') |
---|
| 1355 | i += 1 |
---|
| 1356 | |
---|
| 1357 | assert lines[i].strip().upper().startswith('CROSS') |
---|
| 1358 | i += 1 |
---|
| 1359 | |
---|
| 1360 | assert lines[i].strip().upper().startswith('MAP PROJECTION:') |
---|
| 1361 | projection = lines[i].strip().split(':')[1] |
---|
| 1362 | i += 1 |
---|
| 1363 | |
---|
| 1364 | assert lines[i].strip().upper().startswith('PROJECTION ZONE:') |
---|
| 1365 | zone = int(lines[i].strip().split(':')[1]) |
---|
| 1366 | i += 1 |
---|
| 1367 | |
---|
| 1368 | assert lines[i].strip().upper().startswith('DATUM:') |
---|
| 1369 | datum = lines[i].strip().split(':')[1] |
---|
| 1370 | i += 1 |
---|
| 1371 | |
---|
| 1372 | assert lines[i].strip().upper().startswith('VERTICAL DATUM:') |
---|
| 1373 | i += 1 |
---|
| 1374 | |
---|
| 1375 | assert lines[i].strip().upper().startswith('NUMBER OF REACHES:') |
---|
| 1376 | i += 1 |
---|
| 1377 | |
---|
| 1378 | assert lines[i].strip().upper().startswith('NUMBER OF CROSS-SECTIONS:') |
---|
| 1379 | number_of_cross_sections = int(lines[i].strip().split(':')[1]) |
---|
| 1380 | i += 1 |
---|
| 1381 | |
---|
| 1382 | |
---|
| 1383 | #Now read all points |
---|
| 1384 | points = [] |
---|
| 1385 | elevation = [] |
---|
| 1386 | for j, entries in enumerate(_read_hecras_cross_sections(lines[i:])): |
---|
| 1387 | for k, entry in enumerate(entries): |
---|
| 1388 | points.append(entry[:2]) |
---|
| 1389 | elevation.append(entry[2]) |
---|
| 1390 | |
---|
| 1391 | |
---|
| 1392 | msg = 'Actual #number_of_cross_sections == %d, Reported as %d'\ |
---|
| 1393 | %(j+1, number_of_cross_sections) |
---|
| 1394 | assert j+1 == number_of_cross_sections, msg |
---|
| 1395 | |
---|
| 1396 | #Get output file |
---|
| 1397 | if basename_out == None: |
---|
| 1398 | ptsname = root + '.pts' |
---|
| 1399 | else: |
---|
| 1400 | ptsname = basename_out + '.pts' |
---|
| 1401 | |
---|
| 1402 | #FIXME (DSG-ON): use loadASCII export_points_file |
---|
| 1403 | if verbose: print 'Store to NetCDF file %s' %ptsname |
---|
| 1404 | # NetCDF file definition |
---|
| 1405 | outfile = NetCDFFile(ptsname, 'w') |
---|
| 1406 | |
---|
| 1407 | #Create new file |
---|
| 1408 | outfile.institution = 'Geoscience Australia' |
---|
| 1409 | outfile.description = 'NetCDF pts format for compact and portable ' +\ |
---|
| 1410 | 'storage of spatial point data derived from HEC-RAS' |
---|
| 1411 | |
---|
| 1412 | #Georeferencing |
---|
| 1413 | outfile.zone = zone |
---|
| 1414 | outfile.xllcorner = 0.0 |
---|
| 1415 | outfile.yllcorner = 0.0 |
---|
| 1416 | outfile.false_easting = 500000 #FIXME: Use defaults from e.g. config.py |
---|
| 1417 | outfile.false_northing = 1000000 |
---|
| 1418 | |
---|
| 1419 | outfile.projection = projection |
---|
| 1420 | outfile.datum = datum |
---|
| 1421 | outfile.units = units |
---|
| 1422 | |
---|
| 1423 | |
---|
| 1424 | outfile.createDimension('number_of_points', len(points)) |
---|
| 1425 | outfile.createDimension('number_of_dimensions', 2) #This is 2d data |
---|
| 1426 | |
---|
| 1427 | # variable definitions |
---|
| 1428 | outfile.createVariable('points', Float, ('number_of_points', |
---|
| 1429 | 'number_of_dimensions')) |
---|
| 1430 | outfile.createVariable('elevation', Float, ('number_of_points',)) |
---|
| 1431 | |
---|
| 1432 | # Get handles to the variables |
---|
| 1433 | outfile.variables['points'][:] = points |
---|
| 1434 | outfile.variables['elevation'][:] = elevation |
---|
| 1435 | |
---|
| 1436 | outfile.close() |
---|
| 1437 | |
---|
| 1438 | |
---|
| 1439 | |
---|
| 1440 | def sww2dem(basename_in, basename_out = None, |
---|
| 1441 | quantity = None, |
---|
| 1442 | timestep = None, |
---|
| 1443 | reduction = None, |
---|
| 1444 | cellsize = 10, |
---|
| 1445 | NODATA_value = -9999, |
---|
| 1446 | easting_min = None, |
---|
| 1447 | easting_max = None, |
---|
| 1448 | northing_min = None, |
---|
| 1449 | northing_max = None, |
---|
| 1450 | expand_search = False, #To avoid intractable situations (This will be fixed when least_squares gets redesigned) |
---|
| 1451 | verbose = False, |
---|
| 1452 | origin = None, |
---|
| 1453 | datum = 'WGS84', |
---|
| 1454 | format = 'ers'): |
---|
| 1455 | |
---|
| 1456 | """Read SWW file and convert to Digitial Elevation model format (.asc or .ers) |
---|
| 1457 | |
---|
| 1458 | Example (ASC): |
---|
| 1459 | |
---|
| 1460 | ncols 3121 |
---|
| 1461 | nrows 1800 |
---|
| 1462 | xllcorner 722000 |
---|
| 1463 | yllcorner 5893000 |
---|
| 1464 | cellsize 25 |
---|
| 1465 | NODATA_value -9999 |
---|
| 1466 | 138.3698 137.4194 136.5062 135.5558 .......... |
---|
| 1467 | |
---|
| 1468 | Also write accompanying file with same basename_in but extension .prj |
---|
| 1469 | used to fix the UTM zone, datum, false northings and eastings. |
---|
| 1470 | |
---|
| 1471 | The prj format is assumed to be as |
---|
| 1472 | |
---|
| 1473 | Projection UTM |
---|
| 1474 | Zone 56 |
---|
| 1475 | Datum WGS84 |
---|
| 1476 | Zunits NO |
---|
| 1477 | Units METERS |
---|
| 1478 | Spheroid WGS84 |
---|
| 1479 | Xshift 0.0000000000 |
---|
| 1480 | Yshift 10000000.0000000000 |
---|
| 1481 | Parameters |
---|
| 1482 | |
---|
| 1483 | |
---|
| 1484 | The parameter quantity must be the name of an existing quantity or |
---|
| 1485 | an expression involving existing quantities. The default is |
---|
| 1486 | 'elevation'. |
---|
| 1487 | |
---|
| 1488 | if timestep (an index) is given, output quantity at that timestep |
---|
| 1489 | |
---|
| 1490 | if reduction is given use that to reduce quantity over all timesteps. |
---|
| 1491 | |
---|
| 1492 | datum |
---|
| 1493 | |
---|
| 1494 | format can be either 'asc' or 'ers' |
---|
| 1495 | """ |
---|
| 1496 | |
---|
| 1497 | import sys |
---|
| 1498 | from Numeric import array, Float, concatenate, NewAxis, zeros, reshape, sometrue |
---|
| 1499 | from Numeric import array2string |
---|
| 1500 | |
---|
| 1501 | from utilities.polygon import inside_polygon, outside_polygon, separate_points_by_polygon |
---|
| 1502 | from util import apply_expression_to_dictionary |
---|
| 1503 | |
---|
| 1504 | msg = 'Format must be either asc or ers' |
---|
| 1505 | assert format.lower() in ['asc', 'ers'], msg |
---|
| 1506 | |
---|
| 1507 | |
---|
| 1508 | false_easting = 500000 |
---|
| 1509 | false_northing = 10000000 |
---|
| 1510 | |
---|
| 1511 | if quantity is None: |
---|
| 1512 | quantity = 'elevation' |
---|
| 1513 | |
---|
| 1514 | if reduction is None: |
---|
| 1515 | reduction = max |
---|
| 1516 | |
---|
| 1517 | if basename_out is None: |
---|
| 1518 | basename_out = basename_in + '_%s' %quantity |
---|
| 1519 | |
---|
| 1520 | swwfile = basename_in + '.sww' |
---|
| 1521 | demfile = basename_out + '.' + format |
---|
| 1522 | # Note the use of a .ers extension is optional (write_ermapper_grid will |
---|
| 1523 | # deal with either option |
---|
| 1524 | |
---|
| 1525 | #Read sww file |
---|
| 1526 | if verbose: print 'Reading from %s' %swwfile |
---|
| 1527 | from Scientific.IO.NetCDF import NetCDFFile |
---|
| 1528 | fid = NetCDFFile(swwfile) |
---|
| 1529 | |
---|
| 1530 | #Get extent and reference |
---|
| 1531 | x = fid.variables['x'][:] |
---|
| 1532 | y = fid.variables['y'][:] |
---|
| 1533 | volumes = fid.variables['volumes'][:] |
---|
| 1534 | |
---|
| 1535 | number_of_timesteps = fid.dimensions['number_of_timesteps'] |
---|
| 1536 | number_of_points = fid.dimensions['number_of_points'] |
---|
| 1537 | if origin is None: |
---|
| 1538 | |
---|
| 1539 | #Get geo_reference |
---|
| 1540 | #sww files don't have to have a geo_ref |
---|
| 1541 | try: |
---|
| 1542 | geo_reference = Geo_reference(NetCDFObject=fid) |
---|
| 1543 | except AttributeError, e: |
---|
| 1544 | geo_reference = Geo_reference() #Default georef object |
---|
| 1545 | |
---|
| 1546 | xllcorner = geo_reference.get_xllcorner() |
---|
| 1547 | yllcorner = geo_reference.get_yllcorner() |
---|
| 1548 | zone = geo_reference.get_zone() |
---|
| 1549 | else: |
---|
| 1550 | zone = origin[0] |
---|
| 1551 | xllcorner = origin[1] |
---|
| 1552 | yllcorner = origin[2] |
---|
| 1553 | |
---|
| 1554 | |
---|
| 1555 | |
---|
| 1556 | #FIXME: Refactor using code from file_function.statistics |
---|
| 1557 | #Something like print swwstats(swwname) |
---|
| 1558 | if verbose: |
---|
| 1559 | x = fid.variables['x'][:] |
---|
| 1560 | y = fid.variables['y'][:] |
---|
| 1561 | times = fid.variables['time'][:] |
---|
| 1562 | print '------------------------------------------------' |
---|
| 1563 | print 'Statistics of SWW file:' |
---|
| 1564 | print ' Name: %s' %swwfile |
---|
| 1565 | print ' Reference:' |
---|
| 1566 | print ' Lower left corner: [%f, %f]'\ |
---|
| 1567 | %(xllcorner, yllcorner) |
---|
| 1568 | print ' Start time: %f' %fid.starttime[0] |
---|
| 1569 | print ' Extent:' |
---|
| 1570 | print ' x [m] in [%f, %f], len(x) == %d'\ |
---|
| 1571 | %(min(x.flat), max(x.flat), len(x.flat)) |
---|
| 1572 | print ' y [m] in [%f, %f], len(y) == %d'\ |
---|
| 1573 | %(min(y.flat), max(y.flat), len(y.flat)) |
---|
| 1574 | print ' t [s] in [%f, %f], len(t) == %d'\ |
---|
| 1575 | %(min(times), max(times), len(times)) |
---|
| 1576 | print ' Quantities [SI units]:' |
---|
| 1577 | for name in ['stage', 'xmomentum', 'ymomentum', 'elevation']: |
---|
| 1578 | q = fid.variables[name][:].flat |
---|
| 1579 | print ' %s in [%f, %f]' %(name, min(q), max(q)) |
---|
| 1580 | |
---|
| 1581 | |
---|
| 1582 | |
---|
| 1583 | |
---|
| 1584 | |
---|
| 1585 | #Get quantity and reduce if applicable |
---|
| 1586 | if verbose: print 'Processing quantity %s' %quantity |
---|
| 1587 | |
---|
| 1588 | #Turn NetCDF objects into Numeric arrays |
---|
| 1589 | quantity_dict = {} |
---|
| 1590 | for name in fid.variables.keys(): |
---|
| 1591 | quantity_dict[name] = fid.variables[name][:] |
---|
| 1592 | |
---|
| 1593 | |
---|
| 1594 | q = apply_expression_to_dictionary(quantity, quantity_dict) |
---|
| 1595 | |
---|
| 1596 | |
---|
| 1597 | |
---|
| 1598 | if len(q.shape) == 2: |
---|
| 1599 | #q has a time component and needs to be reduced along |
---|
| 1600 | #the temporal dimension |
---|
| 1601 | if verbose: print 'Reducing quantity %s' %quantity |
---|
| 1602 | q_reduced = zeros( number_of_points, Float ) |
---|
| 1603 | |
---|
| 1604 | for k in range(number_of_points): |
---|
| 1605 | q_reduced[k] = reduction( q[:,k] ) |
---|
| 1606 | |
---|
| 1607 | q = q_reduced |
---|
| 1608 | |
---|
| 1609 | #Post condition: Now q has dimension: number_of_points |
---|
| 1610 | assert len(q.shape) == 1 |
---|
| 1611 | assert q.shape[0] == number_of_points |
---|
| 1612 | |
---|
| 1613 | |
---|
| 1614 | if verbose: |
---|
| 1615 | print 'Processed values for %s are in [%f, %f]' %(quantity, min(q), max(q)) |
---|
| 1616 | |
---|
| 1617 | |
---|
| 1618 | #Create grid and update xll/yll corner and x,y |
---|
| 1619 | |
---|
| 1620 | #Relative extent |
---|
| 1621 | if easting_min is None: |
---|
| 1622 | xmin = min(x) |
---|
| 1623 | else: |
---|
| 1624 | xmin = easting_min - xllcorner |
---|
| 1625 | |
---|
| 1626 | if easting_max is None: |
---|
| 1627 | xmax = max(x) |
---|
| 1628 | else: |
---|
| 1629 | xmax = easting_max - xllcorner |
---|
| 1630 | |
---|
| 1631 | if northing_min is None: |
---|
| 1632 | ymin = min(y) |
---|
| 1633 | else: |
---|
| 1634 | ymin = northing_min - yllcorner |
---|
| 1635 | |
---|
| 1636 | if northing_max is None: |
---|
| 1637 | ymax = max(y) |
---|
| 1638 | else: |
---|
| 1639 | ymax = northing_max - yllcorner |
---|
| 1640 | |
---|
| 1641 | |
---|
| 1642 | |
---|
| 1643 | if verbose: print 'Creating grid' |
---|
| 1644 | ncols = int((xmax-xmin)/cellsize)+1 |
---|
| 1645 | nrows = int((ymax-ymin)/cellsize)+1 |
---|
| 1646 | |
---|
| 1647 | |
---|
| 1648 | #New absolute reference and coordinates |
---|
| 1649 | newxllcorner = xmin+xllcorner |
---|
| 1650 | newyllcorner = ymin+yllcorner |
---|
| 1651 | |
---|
| 1652 | x = x+xllcorner-newxllcorner |
---|
| 1653 | y = y+yllcorner-newyllcorner |
---|
| 1654 | |
---|
| 1655 | vertex_points = concatenate ((x[:, NewAxis] ,y[:, NewAxis]), axis = 1) |
---|
| 1656 | assert len(vertex_points.shape) == 2 |
---|
| 1657 | |
---|
| 1658 | |
---|
| 1659 | |
---|
| 1660 | grid_points = zeros ( (ncols*nrows, 2), Float ) |
---|
| 1661 | |
---|
| 1662 | |
---|
| 1663 | for i in xrange(nrows): |
---|
| 1664 | if format.lower() == 'asc': |
---|
| 1665 | yg = i*cellsize |
---|
| 1666 | else: |
---|
| 1667 | #this will flip the order of the y values for ers |
---|
| 1668 | yg = (nrows-i)*cellsize |
---|
| 1669 | |
---|
| 1670 | for j in xrange(ncols): |
---|
| 1671 | xg = j*cellsize |
---|
| 1672 | k = i*ncols + j |
---|
| 1673 | |
---|
| 1674 | grid_points[k,0] = xg |
---|
| 1675 | grid_points[k,1] = yg |
---|
| 1676 | |
---|
| 1677 | #Interpolate |
---|
| 1678 | from least_squares import Interpolation |
---|
| 1679 | |
---|
| 1680 | |
---|
| 1681 | #FIXME: This should be done with precrop = True (?), otherwise it'll |
---|
| 1682 | #take forever. With expand_search set to False, some grid points might |
---|
| 1683 | #miss out.... This will be addressed though Duncan's refactoring of least_squares |
---|
| 1684 | |
---|
| 1685 | interp = Interpolation(vertex_points, volumes, grid_points, alpha=0.0, |
---|
| 1686 | precrop = False, |
---|
| 1687 | expand_search = expand_search, |
---|
| 1688 | verbose = verbose) |
---|
| 1689 | |
---|
| 1690 | |
---|
| 1691 | |
---|
| 1692 | #Interpolate using quantity values |
---|
| 1693 | if verbose: print 'Interpolating' |
---|
| 1694 | grid_values = interp.interpolate(q).flat |
---|
| 1695 | |
---|
| 1696 | |
---|
| 1697 | if verbose: |
---|
| 1698 | print 'Interpolated values are in [%f, %f]' %(min(grid_values), |
---|
| 1699 | max(grid_values)) |
---|
| 1700 | |
---|
| 1701 | #Assign NODATA_value to all points outside bounding polygon (from interpolation mesh) |
---|
| 1702 | P = interp.mesh.get_boundary_polygon() |
---|
| 1703 | outside_indices = outside_polygon(grid_points, P, closed=True) |
---|
| 1704 | |
---|
| 1705 | for i in outside_indices: |
---|
| 1706 | grid_values[i] = NODATA_value |
---|
| 1707 | |
---|
| 1708 | |
---|
| 1709 | |
---|
| 1710 | |
---|
| 1711 | if format.lower() == 'ers': |
---|
| 1712 | # setup ERS header information |
---|
| 1713 | grid_values = reshape(grid_values,(nrows, ncols)) |
---|
| 1714 | header = {} |
---|
| 1715 | header['datum'] = '"' + datum + '"' |
---|
| 1716 | # FIXME The use of hardwired UTM and zone number needs to be made optional |
---|
| 1717 | # FIXME Also need an automatic test for coordinate type (i.e. EN or LL) |
---|
| 1718 | header['projection'] = '"UTM-' + str(zone) + '"' |
---|
| 1719 | header['coordinatetype'] = 'EN' |
---|
| 1720 | if header['coordinatetype'] == 'LL': |
---|
| 1721 | header['longitude'] = str(newxllcorner) |
---|
| 1722 | header['latitude'] = str(newyllcorner) |
---|
| 1723 | elif header['coordinatetype'] == 'EN': |
---|
| 1724 | header['eastings'] = str(newxllcorner) |
---|
| 1725 | header['northings'] = str(newyllcorner) |
---|
| 1726 | header['nullcellvalue'] = str(NODATA_value) |
---|
| 1727 | header['xdimension'] = str(cellsize) |
---|
| 1728 | header['ydimension'] = str(cellsize) |
---|
| 1729 | header['value'] = '"' + quantity + '"' |
---|
| 1730 | #header['celltype'] = 'IEEE8ByteReal' #FIXME: Breaks unit test |
---|
| 1731 | |
---|
| 1732 | |
---|
| 1733 | #Write |
---|
| 1734 | if verbose: print 'Writing %s' %demfile |
---|
| 1735 | import ermapper_grids |
---|
| 1736 | ermapper_grids.write_ermapper_grid(demfile, grid_values, header) |
---|
| 1737 | |
---|
| 1738 | fid.close() |
---|
| 1739 | else: |
---|
| 1740 | #Write to Ascii format |
---|
| 1741 | |
---|
| 1742 | #Write prj file |
---|
| 1743 | prjfile = basename_out + '.prj' |
---|
| 1744 | |
---|
| 1745 | if verbose: print 'Writing %s' %prjfile |
---|
| 1746 | prjid = open(prjfile, 'w') |
---|
| 1747 | prjid.write('Projection %s\n' %'UTM') |
---|
| 1748 | prjid.write('Zone %d\n' %zone) |
---|
| 1749 | prjid.write('Datum %s\n' %datum) |
---|
| 1750 | prjid.write('Zunits NO\n') |
---|
| 1751 | prjid.write('Units METERS\n') |
---|
| 1752 | prjid.write('Spheroid %s\n' %datum) |
---|
| 1753 | prjid.write('Xshift %d\n' %false_easting) |
---|
| 1754 | prjid.write('Yshift %d\n' %false_northing) |
---|
| 1755 | prjid.write('Parameters\n') |
---|
| 1756 | prjid.close() |
---|
| 1757 | |
---|
| 1758 | |
---|
| 1759 | |
---|
| 1760 | if verbose: print 'Writing %s' %demfile |
---|
| 1761 | |
---|
| 1762 | ascid = open(demfile, 'w') |
---|
| 1763 | |
---|
| 1764 | ascid.write('ncols %d\n' %ncols) |
---|
| 1765 | ascid.write('nrows %d\n' %nrows) |
---|
| 1766 | ascid.write('xllcorner %d\n' %newxllcorner) |
---|
| 1767 | ascid.write('yllcorner %d\n' %newyllcorner) |
---|
| 1768 | ascid.write('cellsize %f\n' %cellsize) |
---|
| 1769 | ascid.write('NODATA_value %d\n' %NODATA_value) |
---|
| 1770 | |
---|
| 1771 | |
---|
| 1772 | #Get bounding polygon from mesh |
---|
| 1773 | #P = interp.mesh.get_boundary_polygon() |
---|
| 1774 | #inside_indices = inside_polygon(grid_points, P) |
---|
| 1775 | |
---|
| 1776 | for i in range(nrows): |
---|
| 1777 | if verbose and i%((nrows+10)/10)==0: |
---|
| 1778 | print 'Doing row %d of %d' %(i, nrows) |
---|
| 1779 | |
---|
| 1780 | base_index = (nrows-i-1)*ncols |
---|
| 1781 | |
---|
| 1782 | slice = grid_values[base_index:base_index+ncols] |
---|
| 1783 | s = array2string(slice, max_line_width=sys.maxint) |
---|
| 1784 | ascid.write(s[1:-1] + '\n') |
---|
| 1785 | |
---|
| 1786 | |
---|
| 1787 | #print |
---|
| 1788 | #for j in range(ncols): |
---|
| 1789 | # index = base_index+j# |
---|
| 1790 | # print grid_values[index], |
---|
| 1791 | # ascid.write('%f ' %grid_values[index]) |
---|
| 1792 | #ascid.write('\n') |
---|
| 1793 | |
---|
| 1794 | #Close |
---|
| 1795 | ascid.close() |
---|
| 1796 | fid.close() |
---|
| 1797 | |
---|
| 1798 | #Backwards compatibility |
---|
| 1799 | def sww2asc(basename_in, basename_out = None, |
---|
| 1800 | quantity = None, |
---|
| 1801 | timestep = None, |
---|
| 1802 | reduction = None, |
---|
| 1803 | cellsize = 10, |
---|
| 1804 | verbose = False, |
---|
| 1805 | origin = None): |
---|
| 1806 | print 'sww2asc will soon be obsoleted - please use sww2dem' |
---|
| 1807 | sww2dem(basename_in, |
---|
| 1808 | basename_out = basename_out, |
---|
| 1809 | quantity = quantity, |
---|
| 1810 | timestep = timestep, |
---|
| 1811 | reduction = reduction, |
---|
| 1812 | cellsize = cellsize, |
---|
| 1813 | verbose = verbose, |
---|
| 1814 | origin = origin, |
---|
| 1815 | datum = 'WGS84', |
---|
| 1816 | format = 'asc') |
---|
| 1817 | |
---|
| 1818 | def sww2ers(basename_in, basename_out = None, |
---|
| 1819 | quantity = None, |
---|
| 1820 | timestep = None, |
---|
| 1821 | reduction = None, |
---|
| 1822 | cellsize = 10, |
---|
| 1823 | verbose = False, |
---|
| 1824 | origin = None, |
---|
| 1825 | datum = 'WGS84'): |
---|
| 1826 | print 'sww2ers will soon be obsoleted - please use sww2dem' |
---|
| 1827 | sww2dem(basename_in, |
---|
| 1828 | basename_out = basename_out, |
---|
| 1829 | quantity = quantity, |
---|
| 1830 | timestep = timestep, |
---|
| 1831 | reduction = reduction, |
---|
| 1832 | cellsize = cellsize, |
---|
| 1833 | verbose = verbose, |
---|
| 1834 | origin = origin, |
---|
| 1835 | datum = datum, |
---|
| 1836 | format = 'ers') |
---|
| 1837 | ################################# END COMPATIBILITY ############## |
---|
| 1838 | |
---|
| 1839 | |
---|
| 1840 | |
---|
| 1841 | |
---|
| 1842 | def convert_dem_from_ascii2netcdf(basename_in, basename_out = None, |
---|
| 1843 | verbose=False): |
---|
| 1844 | """Read Digitial Elevation model from the following ASCII format (.asc) |
---|
| 1845 | |
---|
| 1846 | Example: |
---|
| 1847 | |
---|
| 1848 | ncols 3121 |
---|
| 1849 | nrows 1800 |
---|
| 1850 | xllcorner 722000 |
---|
| 1851 | yllcorner 5893000 |
---|
| 1852 | cellsize 25 |
---|
| 1853 | NODATA_value -9999 |
---|
| 1854 | 138.3698 137.4194 136.5062 135.5558 .......... |
---|
| 1855 | |
---|
| 1856 | Convert basename_in + '.asc' to NetCDF format (.dem) |
---|
| 1857 | mimicking the ASCII format closely. |
---|
| 1858 | |
---|
| 1859 | |
---|
| 1860 | An accompanying file with same basename_in but extension .prj must exist |
---|
| 1861 | and is used to fix the UTM zone, datum, false northings and eastings. |
---|
| 1862 | |
---|
| 1863 | The prj format is assumed to be as |
---|
| 1864 | |
---|
| 1865 | Projection UTM |
---|
| 1866 | Zone 56 |
---|
| 1867 | Datum WGS84 |
---|
| 1868 | Zunits NO |
---|
| 1869 | Units METERS |
---|
| 1870 | Spheroid WGS84 |
---|
| 1871 | Xshift 0.0000000000 |
---|
| 1872 | Yshift 10000000.0000000000 |
---|
| 1873 | Parameters |
---|
| 1874 | """ |
---|
| 1875 | |
---|
| 1876 | import os |
---|
| 1877 | from Scientific.IO.NetCDF import NetCDFFile |
---|
| 1878 | from Numeric import Float, array |
---|
| 1879 | |
---|
| 1880 | #root, ext = os.path.splitext(basename_in) |
---|
| 1881 | root = basename_in |
---|
| 1882 | |
---|
| 1883 | ########################################### |
---|
| 1884 | # Read Meta data |
---|
| 1885 | if verbose: print 'Reading METADATA from %s' %root + '.prj' |
---|
| 1886 | metadatafile = open(root + '.prj') |
---|
| 1887 | metalines = metadatafile.readlines() |
---|
| 1888 | metadatafile.close() |
---|
| 1889 | |
---|
| 1890 | L = metalines[0].strip().split() |
---|
| 1891 | assert L[0].strip().lower() == 'projection' |
---|
| 1892 | projection = L[1].strip() #TEXT |
---|
| 1893 | |
---|
| 1894 | L = metalines[1].strip().split() |
---|
| 1895 | assert L[0].strip().lower() == 'zone' |
---|
| 1896 | zone = int(L[1].strip()) |
---|
| 1897 | |
---|
| 1898 | L = metalines[2].strip().split() |
---|
| 1899 | assert L[0].strip().lower() == 'datum' |
---|
| 1900 | datum = L[1].strip() #TEXT |
---|
| 1901 | |
---|
| 1902 | L = metalines[3].strip().split() |
---|
| 1903 | assert L[0].strip().lower() == 'zunits' #IGNORE |
---|
| 1904 | zunits = L[1].strip() #TEXT |
---|
| 1905 | |
---|
| 1906 | L = metalines[4].strip().split() |
---|
| 1907 | assert L[0].strip().lower() == 'units' |
---|
| 1908 | units = L[1].strip() #TEXT |
---|
| 1909 | |
---|
| 1910 | L = metalines[5].strip().split() |
---|
| 1911 | assert L[0].strip().lower() == 'spheroid' #IGNORE |
---|
| 1912 | spheroid = L[1].strip() #TEXT |
---|
| 1913 | |
---|
| 1914 | L = metalines[6].strip().split() |
---|
| 1915 | assert L[0].strip().lower() == 'xshift' |
---|
| 1916 | false_easting = float(L[1].strip()) |
---|
| 1917 | |
---|
| 1918 | L = metalines[7].strip().split() |
---|
| 1919 | assert L[0].strip().lower() == 'yshift' |
---|
| 1920 | false_northing = float(L[1].strip()) |
---|
| 1921 | |
---|
| 1922 | #print false_easting, false_northing, zone, datum |
---|
| 1923 | |
---|
| 1924 | |
---|
| 1925 | ########################################### |
---|
| 1926 | #Read DEM data |
---|
| 1927 | |
---|
| 1928 | datafile = open(basename_in + '.asc') |
---|
| 1929 | |
---|
| 1930 | if verbose: print 'Reading DEM from %s' %(basename_in + '.asc') |
---|
| 1931 | lines = datafile.readlines() |
---|
| 1932 | datafile.close() |
---|
| 1933 | |
---|
| 1934 | if verbose: print 'Got', len(lines), ' lines' |
---|
| 1935 | |
---|
| 1936 | ncols = int(lines[0].split()[1].strip()) |
---|
| 1937 | nrows = int(lines[1].split()[1].strip()) |
---|
| 1938 | xllcorner = float(lines[2].split()[1].strip()) |
---|
| 1939 | yllcorner = float(lines[3].split()[1].strip()) |
---|
| 1940 | cellsize = float(lines[4].split()[1].strip()) |
---|
| 1941 | NODATA_value = int(lines[5].split()[1].strip()) |
---|
| 1942 | |
---|
| 1943 | assert len(lines) == nrows + 6 |
---|
| 1944 | |
---|
| 1945 | |
---|
| 1946 | ########################################## |
---|
| 1947 | |
---|
| 1948 | |
---|
| 1949 | if basename_out == None: |
---|
| 1950 | netcdfname = root + '.dem' |
---|
| 1951 | else: |
---|
| 1952 | netcdfname = basename_out + '.dem' |
---|
| 1953 | |
---|
| 1954 | if verbose: print 'Store to NetCDF file %s' %netcdfname |
---|
| 1955 | # NetCDF file definition |
---|
| 1956 | fid = NetCDFFile(netcdfname, 'w') |
---|
| 1957 | |
---|
| 1958 | #Create new file |
---|
| 1959 | fid.institution = 'Geoscience Australia' |
---|
| 1960 | fid.description = 'NetCDF DEM format for compact and portable storage ' +\ |
---|
| 1961 | 'of spatial point data' |
---|
| 1962 | |
---|
| 1963 | fid.ncols = ncols |
---|
| 1964 | fid.nrows = nrows |
---|
| 1965 | fid.xllcorner = xllcorner |
---|
| 1966 | fid.yllcorner = yllcorner |
---|
| 1967 | fid.cellsize = cellsize |
---|
| 1968 | fid.NODATA_value = NODATA_value |
---|
| 1969 | |
---|
| 1970 | fid.zone = zone |
---|
| 1971 | fid.false_easting = false_easting |
---|
| 1972 | fid.false_northing = false_northing |
---|
| 1973 | fid.projection = projection |
---|
| 1974 | fid.datum = datum |
---|
| 1975 | fid.units = units |
---|
| 1976 | |
---|
| 1977 | |
---|
| 1978 | # dimension definitions |
---|
| 1979 | fid.createDimension('number_of_rows', nrows) |
---|
| 1980 | fid.createDimension('number_of_columns', ncols) |
---|
| 1981 | |
---|
| 1982 | # variable definitions |
---|
| 1983 | fid.createVariable('elevation', Float, ('number_of_rows', |
---|
| 1984 | 'number_of_columns')) |
---|
| 1985 | |
---|
| 1986 | # Get handles to the variables |
---|
| 1987 | elevation = fid.variables['elevation'] |
---|
| 1988 | |
---|
| 1989 | #Store data |
---|
| 1990 | n = len(lines[6:]) |
---|
| 1991 | for i, line in enumerate(lines[6:]): |
---|
| 1992 | fields = line.split() |
---|
| 1993 | if verbose and i%((n+10)/10)==0: |
---|
| 1994 | print 'Processing row %d of %d' %(i, nrows) |
---|
| 1995 | |
---|
| 1996 | elevation[i, :] = array([float(x) for x in fields]) |
---|
| 1997 | |
---|
| 1998 | fid.close() |
---|
| 1999 | |
---|
| 2000 | |
---|
| 2001 | |
---|
| 2002 | |
---|
| 2003 | |
---|
| 2004 | def ferret2sww(basename_in, basename_out = None, |
---|
| 2005 | verbose = False, |
---|
| 2006 | minlat = None, maxlat = None, |
---|
| 2007 | minlon = None, maxlon = None, |
---|
| 2008 | mint = None, maxt = None, mean_stage = 0, |
---|
| 2009 | origin = None, zscale = 1, |
---|
| 2010 | fail_on_NaN = True, |
---|
| 2011 | NaN_filler = 0, |
---|
| 2012 | elevation = None, |
---|
| 2013 | inverted_bathymetry = False |
---|
| 2014 | ): #FIXME: Bathymetry should be obtained |
---|
| 2015 | #from MOST somehow. |
---|
| 2016 | #Alternatively from elsewhere |
---|
| 2017 | #or, as a last resort, |
---|
| 2018 | #specified here. |
---|
| 2019 | #The value of -100 will work |
---|
| 2020 | #for the Wollongong tsunami |
---|
| 2021 | #scenario but is very hacky |
---|
| 2022 | """Convert 'Ferret' NetCDF format for wave propagation to |
---|
| 2023 | sww format native to pyvolution. |
---|
| 2024 | |
---|
| 2025 | Specify only basename_in and read files of the form |
---|
| 2026 | basefilename_ha.nc, basefilename_ua.nc, basefilename_va.nc containing |
---|
| 2027 | relative height, x-velocity and y-velocity, respectively. |
---|
| 2028 | |
---|
| 2029 | Also convert latitude and longitude to UTM. All coordinates are |
---|
| 2030 | assumed to be given in the GDA94 datum. |
---|
| 2031 | |
---|
| 2032 | min's and max's: If omitted - full extend is used. |
---|
| 2033 | To include a value min may equal it, while max must exceed it. |
---|
| 2034 | Lat and lon are assuemd to be in decimal degrees |
---|
| 2035 | |
---|
| 2036 | origin is a 3-tuple with geo referenced |
---|
| 2037 | UTM coordinates (zone, easting, northing) |
---|
| 2038 | |
---|
| 2039 | nc format has values organised as HA[TIME, LATITUDE, LONGITUDE] |
---|
| 2040 | which means that longitude is the fastest |
---|
| 2041 | varying dimension (row major order, so to speak) |
---|
| 2042 | |
---|
| 2043 | ferret2sww uses grid points as vertices in a triangular grid |
---|
| 2044 | counting vertices from lower left corner upwards, then right |
---|
| 2045 | """ |
---|
| 2046 | |
---|
| 2047 | import os |
---|
| 2048 | from Scientific.IO.NetCDF import NetCDFFile |
---|
| 2049 | from Numeric import Float, Int, Int32, searchsorted, zeros, array |
---|
| 2050 | from Numeric import allclose, around |
---|
| 2051 | |
---|
| 2052 | precision = Float |
---|
| 2053 | |
---|
| 2054 | |
---|
| 2055 | #Get NetCDF data |
---|
| 2056 | if verbose: print 'Reading files %s_*.nc' %basename_in |
---|
| 2057 | file_h = NetCDFFile(basename_in + '_ha.nc', 'r') #Wave amplitude (cm) |
---|
| 2058 | file_u = NetCDFFile(basename_in + '_ua.nc', 'r') #Velocity (x) (cm/s) |
---|
| 2059 | file_v = NetCDFFile(basename_in + '_va.nc', 'r') #Velocity (y) (cm/s) |
---|
| 2060 | file_e = NetCDFFile(basename_in + '_e.nc', 'r') #Elevation (z) (m) |
---|
| 2061 | |
---|
| 2062 | if basename_out is None: |
---|
| 2063 | swwname = basename_in + '.sww' |
---|
| 2064 | else: |
---|
| 2065 | swwname = basename_out + '.sww' |
---|
| 2066 | |
---|
| 2067 | times = file_h.variables['TIME'] |
---|
| 2068 | latitudes = file_h.variables['LAT'] |
---|
| 2069 | longitudes = file_h.variables['LON'] |
---|
| 2070 | |
---|
| 2071 | |
---|
| 2072 | |
---|
| 2073 | #Precision used by most for lat/lon is 4 or 5 decimals |
---|
| 2074 | e_lat = around(file_e.variables['LAT'][:], 5) |
---|
| 2075 | e_lon = around(file_e.variables['LON'][:], 5) |
---|
| 2076 | |
---|
| 2077 | #Check that files are compatible |
---|
| 2078 | assert allclose(latitudes, file_u.variables['LAT']) |
---|
| 2079 | assert allclose(latitudes, file_v.variables['LAT']) |
---|
| 2080 | assert allclose(latitudes, e_lat) |
---|
| 2081 | |
---|
| 2082 | assert allclose(longitudes, file_u.variables['LON']) |
---|
| 2083 | assert allclose(longitudes, file_v.variables['LON']) |
---|
| 2084 | assert allclose(longitudes, e_lon) |
---|
| 2085 | |
---|
| 2086 | |
---|
| 2087 | |
---|
| 2088 | if mint == None: |
---|
| 2089 | jmin = 0 |
---|
| 2090 | else: |
---|
| 2091 | jmin = searchsorted(times, mint) |
---|
| 2092 | |
---|
| 2093 | if maxt == None: |
---|
| 2094 | jmax=len(times) |
---|
| 2095 | else: |
---|
| 2096 | jmax = searchsorted(times, maxt) |
---|
| 2097 | |
---|
| 2098 | if minlat == None: |
---|
| 2099 | kmin=0 |
---|
| 2100 | else: |
---|
| 2101 | kmin = searchsorted(latitudes, minlat) |
---|
| 2102 | |
---|
| 2103 | if maxlat == None: |
---|
| 2104 | kmax = len(latitudes) |
---|
| 2105 | else: |
---|
| 2106 | kmax = searchsorted(latitudes, maxlat) |
---|
| 2107 | |
---|
| 2108 | if minlon == None: |
---|
| 2109 | lmin=0 |
---|
| 2110 | else: |
---|
| 2111 | lmin = searchsorted(longitudes, minlon) |
---|
| 2112 | |
---|
| 2113 | if maxlon == None: |
---|
| 2114 | lmax = len(longitudes) |
---|
| 2115 | else: |
---|
| 2116 | lmax = searchsorted(longitudes, maxlon) |
---|
| 2117 | |
---|
| 2118 | |
---|
| 2119 | |
---|
| 2120 | times = times[jmin:jmax] |
---|
| 2121 | latitudes = latitudes[kmin:kmax] |
---|
| 2122 | longitudes = longitudes[lmin:lmax] |
---|
| 2123 | |
---|
| 2124 | |
---|
| 2125 | if verbose: print 'cropping' |
---|
| 2126 | zname = 'ELEVATION' |
---|
| 2127 | |
---|
| 2128 | |
---|
| 2129 | amplitudes = file_h.variables['HA'][jmin:jmax, kmin:kmax, lmin:lmax] |
---|
| 2130 | uspeed = file_u.variables['UA'][jmin:jmax, kmin:kmax, lmin:lmax] #Lon |
---|
| 2131 | vspeed = file_v.variables['VA'][jmin:jmax, kmin:kmax, lmin:lmax] #Lat |
---|
| 2132 | elevations = file_e.variables[zname][kmin:kmax, lmin:lmax] |
---|
| 2133 | |
---|
| 2134 | # if latitudes2[0]==latitudes[0] and latitudes2[-1]==latitudes[-1]: |
---|
| 2135 | # elevations = file_e.variables['ELEVATION'][kmin:kmax, lmin:lmax] |
---|
| 2136 | # elif latitudes2[0]==latitudes[-1] and latitudes2[-1]==latitudes[0]: |
---|
| 2137 | # from Numeric import asarray |
---|
| 2138 | # elevations=elevations.tolist() |
---|
| 2139 | # elevations.reverse() |
---|
| 2140 | # elevations=asarray(elevations) |
---|
| 2141 | # else: |
---|
| 2142 | # from Numeric import asarray |
---|
| 2143 | # elevations=elevations.tolist() |
---|
| 2144 | # elevations.reverse() |
---|
| 2145 | # elevations=asarray(elevations) |
---|
| 2146 | # 'print hmmm' |
---|
| 2147 | |
---|
| 2148 | |
---|
| 2149 | |
---|
| 2150 | #Get missing values |
---|
| 2151 | nan_ha = file_h.variables['HA'].missing_value[0] |
---|
| 2152 | nan_ua = file_u.variables['UA'].missing_value[0] |
---|
| 2153 | nan_va = file_v.variables['VA'].missing_value[0] |
---|
| 2154 | if hasattr(file_e.variables[zname],'missing_value'): |
---|
| 2155 | nan_e = file_e.variables[zname].missing_value[0] |
---|
| 2156 | else: |
---|
| 2157 | nan_e = None |
---|
| 2158 | |
---|
| 2159 | #Cleanup |
---|
| 2160 | from Numeric import sometrue |
---|
| 2161 | |
---|
| 2162 | missing = (amplitudes == nan_ha) |
---|
| 2163 | if sometrue (missing): |
---|
| 2164 | if fail_on_NaN: |
---|
| 2165 | msg = 'NetCDFFile %s contains missing values'\ |
---|
| 2166 | %(basename_in+'_ha.nc') |
---|
| 2167 | raise msg |
---|
| 2168 | else: |
---|
| 2169 | amplitudes = amplitudes*(missing==0) + missing*NaN_filler |
---|
| 2170 | |
---|
| 2171 | missing = (uspeed == nan_ua) |
---|
| 2172 | if sometrue (missing): |
---|
| 2173 | if fail_on_NaN: |
---|
| 2174 | msg = 'NetCDFFile %s contains missing values'\ |
---|
| 2175 | %(basename_in+'_ua.nc') |
---|
| 2176 | raise msg |
---|
| 2177 | else: |
---|
| 2178 | uspeed = uspeed*(missing==0) + missing*NaN_filler |
---|
| 2179 | |
---|
| 2180 | missing = (vspeed == nan_va) |
---|
| 2181 | if sometrue (missing): |
---|
| 2182 | if fail_on_NaN: |
---|
| 2183 | msg = 'NetCDFFile %s contains missing values'\ |
---|
| 2184 | %(basename_in+'_va.nc') |
---|
| 2185 | raise msg |
---|
| 2186 | else: |
---|
| 2187 | vspeed = vspeed*(missing==0) + missing*NaN_filler |
---|
| 2188 | |
---|
| 2189 | |
---|
| 2190 | missing = (elevations == nan_e) |
---|
| 2191 | if sometrue (missing): |
---|
| 2192 | if fail_on_NaN: |
---|
| 2193 | msg = 'NetCDFFile %s contains missing values'\ |
---|
| 2194 | %(basename_in+'_e.nc') |
---|
| 2195 | raise msg |
---|
| 2196 | else: |
---|
| 2197 | elevations = elevations*(missing==0) + missing*NaN_filler |
---|
| 2198 | |
---|
| 2199 | ####### |
---|
| 2200 | |
---|
| 2201 | |
---|
| 2202 | |
---|
| 2203 | number_of_times = times.shape[0] |
---|
| 2204 | number_of_latitudes = latitudes.shape[0] |
---|
| 2205 | number_of_longitudes = longitudes.shape[0] |
---|
| 2206 | |
---|
| 2207 | assert amplitudes.shape[0] == number_of_times |
---|
| 2208 | assert amplitudes.shape[1] == number_of_latitudes |
---|
| 2209 | assert amplitudes.shape[2] == number_of_longitudes |
---|
| 2210 | |
---|
| 2211 | if verbose: |
---|
| 2212 | print '------------------------------------------------' |
---|
| 2213 | print 'Statistics:' |
---|
| 2214 | print ' Extent (lat/lon):' |
---|
| 2215 | print ' lat in [%f, %f], len(lat) == %d'\ |
---|
| 2216 | %(min(latitudes.flat), max(latitudes.flat), |
---|
| 2217 | len(latitudes.flat)) |
---|
| 2218 | print ' lon in [%f, %f], len(lon) == %d'\ |
---|
| 2219 | %(min(longitudes.flat), max(longitudes.flat), |
---|
| 2220 | len(longitudes.flat)) |
---|
| 2221 | print ' t in [%f, %f], len(t) == %d'\ |
---|
| 2222 | %(min(times.flat), max(times.flat), len(times.flat)) |
---|
| 2223 | |
---|
| 2224 | q = amplitudes.flat |
---|
| 2225 | name = 'Amplitudes (ha) [cm]' |
---|
| 2226 | print ' %s in [%f, %f]' %(name, min(q), max(q)) |
---|
| 2227 | |
---|
| 2228 | q = uspeed.flat |
---|
| 2229 | name = 'Speeds (ua) [cm/s]' |
---|
| 2230 | print ' %s in [%f, %f]' %(name, min(q), max(q)) |
---|
| 2231 | |
---|
| 2232 | q = vspeed.flat |
---|
| 2233 | name = 'Speeds (va) [cm/s]' |
---|
| 2234 | print ' %s in [%f, %f]' %(name, min(q), max(q)) |
---|
| 2235 | |
---|
| 2236 | q = elevations.flat |
---|
| 2237 | name = 'Elevations (e) [m]' |
---|
| 2238 | print ' %s in [%f, %f]' %(name, min(q), max(q)) |
---|
| 2239 | |
---|
| 2240 | |
---|
| 2241 | #print number_of_latitudes, number_of_longitudes |
---|
| 2242 | number_of_points = number_of_latitudes*number_of_longitudes |
---|
| 2243 | number_of_volumes = (number_of_latitudes-1)*(number_of_longitudes-1)*2 |
---|
| 2244 | |
---|
| 2245 | |
---|
| 2246 | file_h.close() |
---|
| 2247 | file_u.close() |
---|
| 2248 | file_v.close() |
---|
| 2249 | file_e.close() |
---|
| 2250 | |
---|
| 2251 | |
---|
| 2252 | # NetCDF file definition |
---|
| 2253 | outfile = NetCDFFile(swwname, 'w') |
---|
| 2254 | |
---|
| 2255 | #Create new file |
---|
| 2256 | outfile.institution = 'Geoscience Australia' |
---|
| 2257 | outfile.description = 'Converted from Ferret files: %s, %s, %s, %s'\ |
---|
| 2258 | %(basename_in + '_ha.nc', |
---|
| 2259 | basename_in + '_ua.nc', |
---|
| 2260 | basename_in + '_va.nc', |
---|
| 2261 | basename_in + '_e.nc') |
---|
| 2262 | |
---|
| 2263 | |
---|
| 2264 | #For sww compatibility |
---|
| 2265 | outfile.smoothing = 'Yes' |
---|
| 2266 | outfile.order = 1 |
---|
| 2267 | |
---|
| 2268 | #Start time in seconds since the epoch (midnight 1/1/1970) |
---|
| 2269 | outfile.starttime = starttime = times[0] |
---|
| 2270 | times = times - starttime #Store relative times |
---|
| 2271 | |
---|
| 2272 | # dimension definitions |
---|
| 2273 | outfile.createDimension('number_of_volumes', number_of_volumes) |
---|
| 2274 | |
---|
| 2275 | outfile.createDimension('number_of_vertices', 3) |
---|
| 2276 | outfile.createDimension('number_of_points', number_of_points) |
---|
| 2277 | |
---|
| 2278 | |
---|
| 2279 | #outfile.createDimension('number_of_timesteps', len(times)) |
---|
| 2280 | outfile.createDimension('number_of_timesteps', len(times)) |
---|
| 2281 | |
---|
| 2282 | # variable definitions |
---|
| 2283 | outfile.createVariable('x', precision, ('number_of_points',)) |
---|
| 2284 | outfile.createVariable('y', precision, ('number_of_points',)) |
---|
| 2285 | outfile.createVariable('elevation', precision, ('number_of_points',)) |
---|
| 2286 | |
---|
| 2287 | #FIXME: Backwards compatibility |
---|
| 2288 | outfile.createVariable('z', precision, ('number_of_points',)) |
---|
| 2289 | ################################# |
---|
| 2290 | |
---|
| 2291 | outfile.createVariable('volumes', Int, ('number_of_volumes', |
---|
| 2292 | 'number_of_vertices')) |
---|
| 2293 | |
---|
| 2294 | outfile.createVariable('time', precision, |
---|
| 2295 | ('number_of_timesteps',)) |
---|
| 2296 | |
---|
| 2297 | outfile.createVariable('stage', precision, |
---|
| 2298 | ('number_of_timesteps', |
---|
| 2299 | 'number_of_points')) |
---|
| 2300 | |
---|
| 2301 | outfile.createVariable('xmomentum', precision, |
---|
| 2302 | ('number_of_timesteps', |
---|
| 2303 | 'number_of_points')) |
---|
| 2304 | |
---|
| 2305 | outfile.createVariable('ymomentum', precision, |
---|
| 2306 | ('number_of_timesteps', |
---|
| 2307 | 'number_of_points')) |
---|
| 2308 | |
---|
| 2309 | |
---|
| 2310 | #Store |
---|
| 2311 | from coordinate_transforms.redfearn import redfearn |
---|
| 2312 | x = zeros(number_of_points, Float) #Easting |
---|
| 2313 | y = zeros(number_of_points, Float) #Northing |
---|
| 2314 | |
---|
| 2315 | |
---|
| 2316 | if verbose: print 'Making triangular grid' |
---|
| 2317 | #Check zone boundaries |
---|
| 2318 | refzone, _, _ = redfearn(latitudes[0],longitudes[0]) |
---|
| 2319 | |
---|
| 2320 | vertices = {} |
---|
| 2321 | i = 0 |
---|
| 2322 | for k, lat in enumerate(latitudes): #Y direction |
---|
| 2323 | for l, lon in enumerate(longitudes): #X direction |
---|
| 2324 | |
---|
| 2325 | vertices[l,k] = i |
---|
| 2326 | |
---|
| 2327 | zone, easting, northing = redfearn(lat,lon) |
---|
| 2328 | |
---|
| 2329 | msg = 'Zone boundary crossed at longitude =', lon |
---|
| 2330 | #assert zone == refzone, msg |
---|
| 2331 | #print '%7.2f %7.2f %8.2f %8.2f' %(lon, lat, easting, northing) |
---|
| 2332 | x[i] = easting |
---|
| 2333 | y[i] = northing |
---|
| 2334 | i += 1 |
---|
| 2335 | |
---|
| 2336 | |
---|
| 2337 | #Construct 2 triangles per 'rectangular' element |
---|
| 2338 | volumes = [] |
---|
| 2339 | for l in range(number_of_longitudes-1): #X direction |
---|
| 2340 | for k in range(number_of_latitudes-1): #Y direction |
---|
| 2341 | v1 = vertices[l,k+1] |
---|
| 2342 | v2 = vertices[l,k] |
---|
| 2343 | v3 = vertices[l+1,k+1] |
---|
| 2344 | v4 = vertices[l+1,k] |
---|
| 2345 | |
---|
| 2346 | volumes.append([v1,v2,v3]) #Upper element |
---|
| 2347 | volumes.append([v4,v3,v2]) #Lower element |
---|
| 2348 | |
---|
| 2349 | volumes = array(volumes) |
---|
| 2350 | |
---|
| 2351 | if origin == None: |
---|
| 2352 | zone = refzone |
---|
| 2353 | xllcorner = min(x) |
---|
| 2354 | yllcorner = min(y) |
---|
| 2355 | else: |
---|
| 2356 | zone = origin[0] |
---|
| 2357 | xllcorner = origin[1] |
---|
| 2358 | yllcorner = origin[2] |
---|
| 2359 | |
---|
| 2360 | |
---|
| 2361 | outfile.xllcorner = xllcorner |
---|
| 2362 | outfile.yllcorner = yllcorner |
---|
| 2363 | outfile.zone = zone |
---|
| 2364 | |
---|
| 2365 | |
---|
| 2366 | if elevation is not None: |
---|
| 2367 | z = elevation |
---|
| 2368 | else: |
---|
| 2369 | if inverted_bathymetry: |
---|
| 2370 | z = -1*elevations |
---|
| 2371 | else: |
---|
| 2372 | z = elevations |
---|
| 2373 | #FIXME: z should be obtained from MOST and passed in here |
---|
| 2374 | |
---|
| 2375 | from Numeric import resize |
---|
| 2376 | z = resize(z,outfile.variables['z'][:].shape) |
---|
| 2377 | outfile.variables['x'][:] = x - xllcorner |
---|
| 2378 | outfile.variables['y'][:] = y - yllcorner |
---|
| 2379 | outfile.variables['z'][:] = z #FIXME HACK |
---|
| 2380 | outfile.variables['elevation'][:] = z |
---|
| 2381 | outfile.variables['time'][:] = times #Store time relative |
---|
| 2382 | outfile.variables['volumes'][:] = volumes.astype(Int32) #On Opteron 64 |
---|
| 2383 | |
---|
| 2384 | |
---|
| 2385 | |
---|
| 2386 | #Time stepping |
---|
| 2387 | stage = outfile.variables['stage'] |
---|
| 2388 | xmomentum = outfile.variables['xmomentum'] |
---|
| 2389 | ymomentum = outfile.variables['ymomentum'] |
---|
| 2390 | |
---|
| 2391 | if verbose: print 'Converting quantities' |
---|
| 2392 | n = len(times) |
---|
| 2393 | for j in range(n): |
---|
| 2394 | if verbose and j%((n+10)/10)==0: print ' Doing %d of %d' %(j, n) |
---|
| 2395 | i = 0 |
---|
| 2396 | for k in range(number_of_latitudes): #Y direction |
---|
| 2397 | for l in range(number_of_longitudes): #X direction |
---|
| 2398 | w = zscale*amplitudes[j,k,l]/100 + mean_stage |
---|
| 2399 | stage[j,i] = w |
---|
| 2400 | h = w - z[i] |
---|
| 2401 | xmomentum[j,i] = uspeed[j,k,l]/100*h |
---|
| 2402 | ymomentum[j,i] = vspeed[j,k,l]/100*h |
---|
| 2403 | i += 1 |
---|
| 2404 | |
---|
| 2405 | #outfile.close() |
---|
| 2406 | |
---|
| 2407 | #FIXME: Refactor using code from file_function.statistics |
---|
| 2408 | #Something like print swwstats(swwname) |
---|
| 2409 | if verbose: |
---|
| 2410 | x = outfile.variables['x'][:] |
---|
| 2411 | y = outfile.variables['y'][:] |
---|
| 2412 | times = outfile.variables['time'][:] |
---|
| 2413 | print '------------------------------------------------' |
---|
| 2414 | print 'Statistics of output file:' |
---|
| 2415 | print ' Name: %s' %swwname |
---|
| 2416 | print ' Reference:' |
---|
| 2417 | print ' Lower left corner: [%f, %f]'\ |
---|
| 2418 | %(xllcorner, yllcorner) |
---|
| 2419 | print ' Start time: %f' %starttime |
---|
| 2420 | print ' Extent:' |
---|
| 2421 | print ' x [m] in [%f, %f], len(x) == %d'\ |
---|
| 2422 | %(min(x.flat), max(x.flat), len(x.flat)) |
---|
| 2423 | print ' y [m] in [%f, %f], len(y) == %d'\ |
---|
| 2424 | %(min(y.flat), max(y.flat), len(y.flat)) |
---|
| 2425 | print ' t [s] in [%f, %f], len(t) == %d'\ |
---|
| 2426 | %(min(times), max(times), len(times)) |
---|
| 2427 | print ' Quantities [SI units]:' |
---|
| 2428 | for name in ['stage', 'xmomentum', 'ymomentum', 'elevation']: |
---|
| 2429 | q = outfile.variables[name][:].flat |
---|
| 2430 | print ' %s in [%f, %f]' %(name, min(q), max(q)) |
---|
| 2431 | |
---|
| 2432 | |
---|
| 2433 | |
---|
| 2434 | outfile.close() |
---|
| 2435 | |
---|
| 2436 | |
---|
| 2437 | |
---|
| 2438 | |
---|
| 2439 | |
---|
| 2440 | def timefile2netcdf(filename, quantity_names = None): |
---|
| 2441 | """Template for converting typical text files with time series to |
---|
| 2442 | NetCDF tms file. |
---|
| 2443 | |
---|
| 2444 | |
---|
| 2445 | The file format is assumed to be either two fields separated by a comma: |
---|
| 2446 | |
---|
| 2447 | time [DD/MM/YY hh:mm:ss], value0 value1 value2 ... |
---|
| 2448 | |
---|
| 2449 | E.g |
---|
| 2450 | |
---|
| 2451 | 31/08/04 00:00:00, 1.328223 0 0 |
---|
| 2452 | 31/08/04 00:15:00, 1.292912 0 0 |
---|
| 2453 | |
---|
| 2454 | will provide a time dependent function f(t) with three attributes |
---|
| 2455 | |
---|
| 2456 | filename is assumed to be the rootname with extenisons .txt and .sww |
---|
| 2457 | """ |
---|
| 2458 | |
---|
| 2459 | import time, calendar |
---|
| 2460 | from Numeric import array |
---|
| 2461 | from config import time_format |
---|
| 2462 | from util import ensure_numeric |
---|
| 2463 | |
---|
| 2464 | |
---|
| 2465 | fid = open(filename + '.txt') |
---|
| 2466 | line = fid.readline() |
---|
| 2467 | fid.close() |
---|
| 2468 | |
---|
| 2469 | fields = line.split(',') |
---|
| 2470 | msg = 'File %s must have the format date, value0 value1 value2 ...' |
---|
| 2471 | assert len(fields) == 2, msg |
---|
| 2472 | |
---|
| 2473 | try: |
---|
| 2474 | starttime = calendar.timegm(time.strptime(fields[0], time_format)) |
---|
| 2475 | except ValueError: |
---|
| 2476 | msg = 'First field in file %s must be' %filename |
---|
| 2477 | msg += ' date-time with format %s.\n' %time_format |
---|
| 2478 | msg += 'I got %s instead.' %fields[0] |
---|
| 2479 | raise msg |
---|
| 2480 | |
---|
| 2481 | |
---|
| 2482 | #Split values |
---|
| 2483 | values = [] |
---|
| 2484 | for value in fields[1].split(): |
---|
| 2485 | values.append(float(value)) |
---|
| 2486 | |
---|
| 2487 | q = ensure_numeric(values) |
---|
| 2488 | |
---|
| 2489 | msg = 'ERROR: File must contain at least one independent value' |
---|
| 2490 | assert len(q.shape) == 1, msg |
---|
| 2491 | |
---|
| 2492 | |
---|
| 2493 | |
---|
| 2494 | #Read times proper |
---|
| 2495 | from Numeric import zeros, Float, alltrue |
---|
| 2496 | from config import time_format |
---|
| 2497 | import time, calendar |
---|
| 2498 | |
---|
| 2499 | fid = open(filename + '.txt') |
---|
| 2500 | lines = fid.readlines() |
---|
| 2501 | fid.close() |
---|
| 2502 | |
---|
| 2503 | N = len(lines) |
---|
| 2504 | d = len(q) |
---|
| 2505 | |
---|
| 2506 | T = zeros(N, Float) #Time |
---|
| 2507 | Q = zeros((N, d), Float) #Values |
---|
| 2508 | |
---|
| 2509 | for i, line in enumerate(lines): |
---|
| 2510 | fields = line.split(',') |
---|
| 2511 | realtime = calendar.timegm(time.strptime(fields[0], time_format)) |
---|
| 2512 | |
---|
| 2513 | T[i] = realtime - starttime |
---|
| 2514 | |
---|
| 2515 | for j, value in enumerate(fields[1].split()): |
---|
| 2516 | Q[i, j] = float(value) |
---|
| 2517 | |
---|
| 2518 | msg = 'File %s must list time as a monotonuosly ' %filename |
---|
| 2519 | msg += 'increasing sequence' |
---|
| 2520 | assert alltrue( T[1:] - T[:-1] > 0 ), msg |
---|
| 2521 | |
---|
| 2522 | |
---|
| 2523 | #Create NetCDF file |
---|
| 2524 | from Scientific.IO.NetCDF import NetCDFFile |
---|
| 2525 | |
---|
| 2526 | fid = NetCDFFile(filename + '.tms', 'w') |
---|
| 2527 | |
---|
| 2528 | |
---|
| 2529 | fid.institution = 'Geoscience Australia' |
---|
| 2530 | fid.description = 'Time series' |
---|
| 2531 | |
---|
| 2532 | |
---|
| 2533 | #Reference point |
---|
| 2534 | #Start time in seconds since the epoch (midnight 1/1/1970) |
---|
| 2535 | #FIXME: Use Georef |
---|
| 2536 | fid.starttime = starttime |
---|
| 2537 | |
---|
| 2538 | |
---|
| 2539 | # dimension definitions |
---|
| 2540 | #fid.createDimension('number_of_volumes', self.number_of_volumes) |
---|
| 2541 | #fid.createDimension('number_of_vertices', 3) |
---|
| 2542 | |
---|
| 2543 | |
---|
| 2544 | fid.createDimension('number_of_timesteps', len(T)) |
---|
| 2545 | |
---|
| 2546 | fid.createVariable('time', Float, ('number_of_timesteps',)) |
---|
| 2547 | |
---|
| 2548 | fid.variables['time'][:] = T |
---|
| 2549 | |
---|
| 2550 | for i in range(Q.shape[1]): |
---|
| 2551 | try: |
---|
| 2552 | name = quantity_names[i] |
---|
| 2553 | except: |
---|
| 2554 | name = 'Attribute%d'%i |
---|
| 2555 | |
---|
| 2556 | fid.createVariable(name, Float, ('number_of_timesteps',)) |
---|
| 2557 | fid.variables[name][:] = Q[:,i] |
---|
| 2558 | |
---|
| 2559 | fid.close() |
---|
| 2560 | |
---|
| 2561 | |
---|
| 2562 | def extent_sww(file_name): |
---|
| 2563 | """ |
---|
| 2564 | Read in an sww file. |
---|
| 2565 | |
---|
| 2566 | Input; |
---|
| 2567 | file_name - the sww file |
---|
| 2568 | |
---|
| 2569 | Output; |
---|
| 2570 | z - Vector of bed elevation |
---|
| 2571 | volumes - Array. Each row has 3 values, representing |
---|
| 2572 | the vertices that define the volume |
---|
| 2573 | time - Vector of the times where there is stage information |
---|
| 2574 | stage - array with respect to time and vertices (x,y) |
---|
| 2575 | """ |
---|
| 2576 | |
---|
| 2577 | |
---|
| 2578 | from Scientific.IO.NetCDF import NetCDFFile |
---|
| 2579 | |
---|
| 2580 | #Check contents |
---|
| 2581 | #Get NetCDF |
---|
| 2582 | fid = NetCDFFile(file_name, 'r') |
---|
| 2583 | |
---|
| 2584 | # Get the variables |
---|
| 2585 | x = fid.variables['x'][:] |
---|
| 2586 | y = fid.variables['y'][:] |
---|
| 2587 | stage = fid.variables['stage'][:] |
---|
| 2588 | #print "stage",stage |
---|
| 2589 | #print "stage.shap",stage.shape |
---|
| 2590 | #print "min(stage.flat), mpythonax(stage.flat)",min(stage.flat), max(stage.flat) |
---|
| 2591 | #print "min(stage)",min(stage) |
---|
| 2592 | |
---|
| 2593 | fid.close() |
---|
| 2594 | |
---|
| 2595 | return [min(x),max(x),min(y),max(y),min(stage.flat),max(stage.flat)] |
---|
| 2596 | |
---|
| 2597 | |
---|
| 2598 | def sww2domain(filename,boundary=None,t=None,\ |
---|
| 2599 | fail_if_NaN=True,NaN_filler=0\ |
---|
| 2600 | ,verbose = False,very_verbose = False): |
---|
| 2601 | """ |
---|
| 2602 | Usage: domain = sww2domain('file.sww',t=time (default = last time in file)) |
---|
| 2603 | |
---|
| 2604 | Boundary is not recommended if domain.smooth is not selected, as it |
---|
| 2605 | uses unique coordinates, but not unique boundaries. This means that |
---|
| 2606 | the boundary file will not be compatable with the coordinates, and will |
---|
| 2607 | give a different final boundary, or crash. |
---|
| 2608 | """ |
---|
| 2609 | NaN=9.969209968386869e+036 |
---|
| 2610 | #initialise NaN. |
---|
| 2611 | |
---|
| 2612 | from Scientific.IO.NetCDF import NetCDFFile |
---|
| 2613 | from shallow_water import Domain |
---|
| 2614 | from Numeric import asarray, transpose, resize |
---|
| 2615 | |
---|
| 2616 | if verbose: print 'Reading from ', filename |
---|
| 2617 | fid = NetCDFFile(filename, 'r') #Open existing file for read |
---|
| 2618 | time = fid.variables['time'] #Timesteps |
---|
| 2619 | if t is None: |
---|
| 2620 | t = time[-1] |
---|
| 2621 | time_interp = get_time_interp(time,t) |
---|
| 2622 | |
---|
| 2623 | # Get the variables as Numeric arrays |
---|
| 2624 | x = fid.variables['x'][:] #x-coordinates of vertices |
---|
| 2625 | y = fid.variables['y'][:] #y-coordinates of vertices |
---|
| 2626 | elevation = fid.variables['elevation'] #Elevation |
---|
| 2627 | stage = fid.variables['stage'] #Water level |
---|
| 2628 | xmomentum = fid.variables['xmomentum'] #Momentum in the x-direction |
---|
| 2629 | ymomentum = fid.variables['ymomentum'] #Momentum in the y-direction |
---|
| 2630 | |
---|
| 2631 | starttime = fid.starttime[0] |
---|
| 2632 | volumes = fid.variables['volumes'][:] #Connectivity |
---|
| 2633 | coordinates=transpose(asarray([x.tolist(),y.tolist()])) |
---|
| 2634 | |
---|
| 2635 | conserved_quantities = [] |
---|
| 2636 | interpolated_quantities = {} |
---|
| 2637 | other_quantities = [] |
---|
| 2638 | |
---|
| 2639 | # get geo_reference |
---|
| 2640 | #sww files don't have to have a geo_ref |
---|
| 2641 | try: |
---|
| 2642 | geo_reference = Geo_reference(NetCDFObject=fid) |
---|
| 2643 | except: #AttributeError, e: |
---|
| 2644 | geo_reference = None |
---|
| 2645 | |
---|
| 2646 | if verbose: print ' getting quantities' |
---|
| 2647 | for quantity in fid.variables.keys(): |
---|
| 2648 | dimensions = fid.variables[quantity].dimensions |
---|
| 2649 | if 'number_of_timesteps' in dimensions: |
---|
| 2650 | conserved_quantities.append(quantity) |
---|
| 2651 | interpolated_quantities[quantity]=\ |
---|
| 2652 | interpolated_quantity(fid.variables[quantity][:],time_interp) |
---|
| 2653 | else: other_quantities.append(quantity) |
---|
| 2654 | |
---|
| 2655 | other_quantities.remove('x') |
---|
| 2656 | other_quantities.remove('y') |
---|
| 2657 | other_quantities.remove('z') |
---|
| 2658 | other_quantities.remove('volumes') |
---|
| 2659 | |
---|
| 2660 | conserved_quantities.remove('time') |
---|
| 2661 | |
---|
| 2662 | if verbose: print ' building domain' |
---|
| 2663 | # From domain.Domain: |
---|
| 2664 | # domain = Domain(coordinates, volumes,\ |
---|
| 2665 | # conserved_quantities = conserved_quantities,\ |
---|
| 2666 | # other_quantities = other_quantities,zone=zone,\ |
---|
| 2667 | # xllcorner=xllcorner, yllcorner=yllcorner) |
---|
| 2668 | |
---|
| 2669 | # From shallow_water.Domain: |
---|
| 2670 | coordinates=coordinates.tolist() |
---|
| 2671 | volumes=volumes.tolist() |
---|
| 2672 | #FIXME:should this be in mesh?(peter row) |
---|
| 2673 | if fid.smoothing == 'Yes': unique = False |
---|
| 2674 | else: unique = True |
---|
| 2675 | if unique: |
---|
| 2676 | coordinates,volumes,boundary=weed(coordinates,volumes,boundary) |
---|
| 2677 | |
---|
| 2678 | |
---|
| 2679 | try: |
---|
| 2680 | domain = Domain(coordinates, volumes, boundary) |
---|
| 2681 | except AssertionError, e: |
---|
| 2682 | fid.close() |
---|
| 2683 | msg = 'Domain could not be created: %s. Perhaps use "fail_if_NaN=False and NaN_filler = ..."' %e |
---|
| 2684 | raise msg |
---|
| 2685 | |
---|
| 2686 | if not boundary is None: |
---|
| 2687 | domain.boundary = boundary |
---|
| 2688 | |
---|
| 2689 | domain.geo_reference = geo_reference |
---|
| 2690 | |
---|
| 2691 | domain.starttime=float(starttime)+float(t) |
---|
| 2692 | domain.time=0.0 |
---|
| 2693 | |
---|
| 2694 | for quantity in other_quantities: |
---|
| 2695 | try: |
---|
| 2696 | NaN = fid.variables[quantity].missing_value |
---|
| 2697 | except: |
---|
| 2698 | pass #quantity has no missing_value number |
---|
| 2699 | X = fid.variables[quantity][:] |
---|
| 2700 | if very_verbose: |
---|
| 2701 | print ' ',quantity |
---|
| 2702 | print ' NaN =',NaN |
---|
| 2703 | print ' max(X)' |
---|
| 2704 | print ' ',max(X) |
---|
| 2705 | print ' max(X)==NaN' |
---|
| 2706 | print ' ',max(X)==NaN |
---|
| 2707 | print '' |
---|
| 2708 | if (max(X)==NaN) or (min(X)==NaN): |
---|
| 2709 | if fail_if_NaN: |
---|
| 2710 | msg = 'quantity "%s" contains no_data entry'%quantity |
---|
| 2711 | raise msg |
---|
| 2712 | else: |
---|
| 2713 | data = (X<>NaN) |
---|
| 2714 | X = (X*data)+(data==0)*NaN_filler |
---|
| 2715 | if unique: |
---|
| 2716 | X = resize(X,(len(X)/3,3)) |
---|
| 2717 | domain.set_quantity(quantity,X) |
---|
| 2718 | # |
---|
| 2719 | for quantity in conserved_quantities: |
---|
| 2720 | try: |
---|
| 2721 | NaN = fid.variables[quantity].missing_value |
---|
| 2722 | except: |
---|
| 2723 | pass #quantity has no missing_value number |
---|
| 2724 | X = interpolated_quantities[quantity] |
---|
| 2725 | if very_verbose: |
---|
| 2726 | print ' ',quantity |
---|
| 2727 | print ' NaN =',NaN |
---|
| 2728 | print ' max(X)' |
---|
| 2729 | print ' ',max(X) |
---|
| 2730 | print ' max(X)==NaN' |
---|
| 2731 | print ' ',max(X)==NaN |
---|
| 2732 | print '' |
---|
| 2733 | if (max(X)==NaN) or (min(X)==NaN): |
---|
| 2734 | if fail_if_NaN: |
---|
| 2735 | msg = 'quantity "%s" contains no_data entry'%quantity |
---|
| 2736 | raise msg |
---|
| 2737 | else: |
---|
| 2738 | data = (X<>NaN) |
---|
| 2739 | X = (X*data)+(data==0)*NaN_filler |
---|
| 2740 | if unique: |
---|
| 2741 | X = resize(X,(X.shape[0]/3,3)) |
---|
| 2742 | domain.set_quantity(quantity,X) |
---|
| 2743 | |
---|
| 2744 | fid.close() |
---|
| 2745 | return domain |
---|
| 2746 | |
---|
| 2747 | def interpolated_quantity(saved_quantity,time_interp): |
---|
| 2748 | |
---|
| 2749 | #given an index and ratio, interpolate quantity with respect to time. |
---|
| 2750 | index,ratio = time_interp |
---|
| 2751 | Q = saved_quantity |
---|
| 2752 | if ratio > 0: |
---|
| 2753 | q = (1-ratio)*Q[index]+ ratio*Q[index+1] |
---|
| 2754 | else: |
---|
| 2755 | q = Q[index] |
---|
| 2756 | #Return vector of interpolated values |
---|
| 2757 | return q |
---|
| 2758 | |
---|
| 2759 | def get_time_interp(time,t=None): |
---|
| 2760 | #Finds the ratio and index for time interpolation. |
---|
| 2761 | #It is borrowed from previous pyvolution code. |
---|
| 2762 | if t is None: |
---|
| 2763 | t=time[-1] |
---|
| 2764 | index = -1 |
---|
| 2765 | ratio = 0. |
---|
| 2766 | else: |
---|
| 2767 | T = time |
---|
| 2768 | tau = t |
---|
| 2769 | index=0 |
---|
| 2770 | msg = 'Time interval derived from file %s [%s:%s]'\ |
---|
| 2771 | %('FIXMEfilename', T[0], T[-1]) |
---|
| 2772 | msg += ' does not match model time: %s' %tau |
---|
| 2773 | if tau < time[0]: raise msg |
---|
| 2774 | if tau > time[-1]: raise msg |
---|
| 2775 | while tau > time[index]: index += 1 |
---|
| 2776 | while tau < time[index]: index -= 1 |
---|
| 2777 | if tau == time[index]: |
---|
| 2778 | #Protect against case where tau == time[-1] (last time) |
---|
| 2779 | # - also works in general when tau == time[i] |
---|
| 2780 | ratio = 0 |
---|
| 2781 | else: |
---|
| 2782 | #t is now between index and index+1 |
---|
| 2783 | ratio = (tau - time[index])/(time[index+1] - time[index]) |
---|
| 2784 | return (index,ratio) |
---|
| 2785 | |
---|
| 2786 | |
---|
| 2787 | def weed(coordinates,volumes,boundary = None): |
---|
| 2788 | if type(coordinates)=='array': |
---|
| 2789 | coordinates = coordinates.tolist() |
---|
| 2790 | if type(volumes)=='array': |
---|
| 2791 | volumes = volumes.tolist() |
---|
| 2792 | |
---|
| 2793 | unique = False |
---|
| 2794 | point_dict = {} |
---|
| 2795 | same_point = {} |
---|
| 2796 | for i in range(len(coordinates)): |
---|
| 2797 | point = tuple(coordinates[i]) |
---|
| 2798 | if point_dict.has_key(point): |
---|
| 2799 | unique = True |
---|
| 2800 | same_point[i]=point |
---|
| 2801 | #to change all point i references to point j |
---|
| 2802 | else: |
---|
| 2803 | point_dict[point]=i |
---|
| 2804 | same_point[i]=point |
---|
| 2805 | |
---|
| 2806 | coordinates = [] |
---|
| 2807 | i = 0 |
---|
| 2808 | for point in point_dict.keys(): |
---|
| 2809 | point = tuple(point) |
---|
| 2810 | coordinates.append(list(point)) |
---|
| 2811 | point_dict[point]=i |
---|
| 2812 | i+=1 |
---|
| 2813 | |
---|
| 2814 | |
---|
| 2815 | for volume in volumes: |
---|
| 2816 | for i in range(len(volume)): |
---|
| 2817 | index = volume[i] |
---|
| 2818 | if index>-1: |
---|
| 2819 | volume[i]=point_dict[same_point[index]] |
---|
| 2820 | |
---|
| 2821 | new_boundary = {} |
---|
| 2822 | if not boundary is None: |
---|
| 2823 | for segment in boundary.keys(): |
---|
| 2824 | point0 = point_dict[same_point[segment[0]]] |
---|
| 2825 | point1 = point_dict[same_point[segment[1]]] |
---|
| 2826 | label = boundary[segment] |
---|
| 2827 | #FIXME should the bounday attributes be concaterated |
---|
| 2828 | #('exterior, pond') or replaced ('pond')(peter row) |
---|
| 2829 | |
---|
| 2830 | if new_boundary.has_key((point0,point1)): |
---|
| 2831 | new_boundary[(point0,point1)]=new_boundary[(point0,point1)]#\ |
---|
| 2832 | #+','+label |
---|
| 2833 | |
---|
| 2834 | elif new_boundary.has_key((point1,point0)): |
---|
| 2835 | new_boundary[(point1,point0)]=new_boundary[(point1,point0)]#\ |
---|
| 2836 | #+','+label |
---|
| 2837 | else: new_boundary[(point0,point1)]=label |
---|
| 2838 | |
---|
| 2839 | boundary = new_boundary |
---|
| 2840 | |
---|
| 2841 | return coordinates,volumes,boundary |
---|
| 2842 | |
---|
| 2843 | |
---|
| 2844 | |
---|
| 2845 | |
---|
| 2846 | |
---|
| 2847 | |
---|
| 2848 | |
---|
| 2849 | |
---|
| 2850 | |
---|
| 2851 | |
---|
| 2852 | ############################################### |
---|
| 2853 | #OBSOLETE STUFF |
---|
| 2854 | #Native checkpoint format. |
---|
| 2855 | #Information needed to recreate a state is preserved |
---|
| 2856 | #FIXME: Rethink and maybe use netcdf format |
---|
| 2857 | def cpt_variable_writer(filename, t, v0, v1, v2): |
---|
| 2858 | """Store all conserved quantities to file |
---|
| 2859 | """ |
---|
| 2860 | |
---|
| 2861 | M, N = v0.shape |
---|
| 2862 | |
---|
| 2863 | FN = create_filename(filename, 'cpt', M, t) |
---|
| 2864 | #print 'Writing to %s' %FN |
---|
| 2865 | |
---|
| 2866 | fid = open(FN, 'w') |
---|
| 2867 | for i in range(M): |
---|
| 2868 | for j in range(N): |
---|
| 2869 | fid.write('%.16e ' %v0[i,j]) |
---|
| 2870 | for j in range(N): |
---|
| 2871 | fid.write('%.16e ' %v1[i,j]) |
---|
| 2872 | for j in range(N): |
---|
| 2873 | fid.write('%.16e ' %v2[i,j]) |
---|
| 2874 | |
---|
| 2875 | fid.write('\n') |
---|
| 2876 | fid.close() |
---|
| 2877 | |
---|
| 2878 | |
---|
| 2879 | def cpt_variable_reader(filename, t, v0, v1, v2): |
---|
| 2880 | """Store all conserved quantities to file |
---|
| 2881 | """ |
---|
| 2882 | |
---|
| 2883 | M, N = v0.shape |
---|
| 2884 | |
---|
| 2885 | FN = create_filename(filename, 'cpt', M, t) |
---|
| 2886 | #print 'Reading from %s' %FN |
---|
| 2887 | |
---|
| 2888 | fid = open(FN) |
---|
| 2889 | |
---|
| 2890 | |
---|
| 2891 | for i in range(M): |
---|
| 2892 | values = fid.readline().split() #Get one line |
---|
| 2893 | |
---|
| 2894 | for j in range(N): |
---|
| 2895 | v0[i,j] = float(values[j]) |
---|
| 2896 | v1[i,j] = float(values[3+j]) |
---|
| 2897 | v2[i,j] = float(values[6+j]) |
---|
| 2898 | |
---|
| 2899 | fid.close() |
---|
| 2900 | |
---|
| 2901 | def cpt_constant_writer(filename, X0, X1, X2, v0, v1, v2): |
---|
| 2902 | """Writes x,y,z,z,z coordinates of triangles constituting the bed |
---|
| 2903 | elevation. |
---|
| 2904 | FIXME: Not in use pt |
---|
| 2905 | """ |
---|
| 2906 | |
---|
| 2907 | M, N = v0.shape |
---|
| 2908 | |
---|
| 2909 | |
---|
| 2910 | print X0 |
---|
| 2911 | import sys; sys.exit() |
---|
| 2912 | FN = create_filename(filename, 'cpt', M) |
---|
| 2913 | print 'Writing to %s' %FN |
---|
| 2914 | |
---|
| 2915 | fid = open(FN, 'w') |
---|
| 2916 | for i in range(M): |
---|
| 2917 | for j in range(2): |
---|
| 2918 | fid.write('%.16e ' %X0[i,j]) #x, y |
---|
| 2919 | for j in range(N): |
---|
| 2920 | fid.write('%.16e ' %v0[i,j]) #z,z,z, |
---|
| 2921 | |
---|
| 2922 | for j in range(2): |
---|
| 2923 | fid.write('%.16e ' %X1[i,j]) #x, y |
---|
| 2924 | for j in range(N): |
---|
| 2925 | fid.write('%.16e ' %v1[i,j]) |
---|
| 2926 | |
---|
| 2927 | for j in range(2): |
---|
| 2928 | fid.write('%.16e ' %X2[i,j]) #x, y |
---|
| 2929 | for j in range(N): |
---|
| 2930 | fid.write('%.16e ' %v2[i,j]) |
---|
| 2931 | |
---|
| 2932 | fid.write('\n') |
---|
| 2933 | fid.close() |
---|
| 2934 | |
---|
| 2935 | |
---|
| 2936 | |
---|
| 2937 | #Function for storing out to e.g. visualisation |
---|
| 2938 | #FIXME: Do we want this? |
---|
| 2939 | #FIXME: Not done yet for this version |
---|
| 2940 | def dat_constant_writer(filename, X0, X1, X2, v0, v1, v2): |
---|
| 2941 | """Writes x,y,z coordinates of triangles constituting the bed elevation. |
---|
| 2942 | """ |
---|
| 2943 | |
---|
| 2944 | M, N = v0.shape |
---|
| 2945 | |
---|
| 2946 | FN = create_filename(filename, 'dat', M) |
---|
| 2947 | #print 'Writing to %s' %FN |
---|
| 2948 | |
---|
| 2949 | fid = open(FN, 'w') |
---|
| 2950 | for i in range(M): |
---|
| 2951 | for j in range(2): |
---|
| 2952 | fid.write('%f ' %X0[i,j]) #x, y |
---|
| 2953 | fid.write('%f ' %v0[i,0]) #z |
---|
| 2954 | |
---|
| 2955 | for j in range(2): |
---|
| 2956 | fid.write('%f ' %X1[i,j]) #x, y |
---|
| 2957 | fid.write('%f ' %v1[i,0]) #z |
---|
| 2958 | |
---|
| 2959 | for j in range(2): |
---|
| 2960 | fid.write('%f ' %X2[i,j]) #x, y |
---|
| 2961 | fid.write('%f ' %v2[i,0]) #z |
---|
| 2962 | |
---|
| 2963 | fid.write('\n') |
---|
| 2964 | fid.close() |
---|
| 2965 | |
---|
| 2966 | |
---|
| 2967 | |
---|
| 2968 | def dat_variable_writer(filename, t, v0, v1, v2): |
---|
| 2969 | """Store water height to file |
---|
| 2970 | """ |
---|
| 2971 | |
---|
| 2972 | M, N = v0.shape |
---|
| 2973 | |
---|
| 2974 | FN = create_filename(filename, 'dat', M, t) |
---|
| 2975 | #print 'Writing to %s' %FN |
---|
| 2976 | |
---|
| 2977 | fid = open(FN, 'w') |
---|
| 2978 | for i in range(M): |
---|
| 2979 | fid.write('%.4f ' %v0[i,0]) |
---|
| 2980 | fid.write('%.4f ' %v1[i,0]) |
---|
| 2981 | fid.write('%.4f ' %v2[i,0]) |
---|
| 2982 | |
---|
| 2983 | fid.write('\n') |
---|
| 2984 | fid.close() |
---|
| 2985 | |
---|
| 2986 | |
---|
| 2987 | def read_sww(filename): |
---|
| 2988 | """Read sww Net CDF file containing Shallow Water Wave simulation |
---|
| 2989 | |
---|
| 2990 | The integer array volumes is of shape Nx3 where N is the number of |
---|
| 2991 | triangles in the mesh. |
---|
| 2992 | |
---|
| 2993 | Each entry in volumes is an index into the x,y arrays (the location). |
---|
| 2994 | |
---|
| 2995 | Quantities stage, elevation, xmomentum and ymomentum are all in arrays of dimensions |
---|
| 2996 | number_of_timesteps, number_of_points. |
---|
| 2997 | |
---|
| 2998 | The momentum is not always stored. |
---|
| 2999 | |
---|
| 3000 | """ |
---|
| 3001 | from Scientific.IO.NetCDF import NetCDFFile |
---|
| 3002 | print 'Reading from ', filename |
---|
| 3003 | fid = NetCDFFile(filename, 'r') #Open existing file for read |
---|
| 3004 | #latitude, longitude |
---|
| 3005 | # Get the variables as Numeric arrays |
---|
| 3006 | x = fid.variables['x'] #x-coordinates of vertices |
---|
| 3007 | y = fid.variables['y'] #y-coordinates of vertices |
---|
| 3008 | z = fid.variables['elevation'] #Elevation |
---|
| 3009 | time = fid.variables['time'] #Timesteps |
---|
| 3010 | stage = fid.variables['stage'] #Water level |
---|
| 3011 | #xmomentum = fid.variables['xmomentum'] #Momentum in the x-direction |
---|
| 3012 | #ymomentum = fid.variables['ymomentum'] #Momentum in the y-direction |
---|
| 3013 | |
---|
| 3014 | volumes = fid.variables['volumes'] #Connectivity |
---|
| 3015 | |
---|
| 3016 | #FIXME (Ole): What is this? |
---|
| 3017 | # Why isn't anything returned? |
---|
| 3018 | # Where's the unit test? |
---|
| 3019 | |
---|
| 3020 | def decimate_dem(basename_in, stencil, cellsize_new, basename_out=None, |
---|
| 3021 | verbose=False): |
---|
| 3022 | """Read Digitial Elevation model from the following NetCDF format (.dem) |
---|
| 3023 | |
---|
| 3024 | Example: |
---|
| 3025 | |
---|
| 3026 | ncols 3121 |
---|
| 3027 | nrows 1800 |
---|
| 3028 | xllcorner 722000 |
---|
| 3029 | yllcorner 5893000 |
---|
| 3030 | cellsize 25 |
---|
| 3031 | NODATA_value -9999 |
---|
| 3032 | 138.3698 137.4194 136.5062 135.5558 .......... |
---|
| 3033 | |
---|
| 3034 | Decimate data to cellsize_new using stencil and write to NetCDF dem format. |
---|
| 3035 | """ |
---|
| 3036 | |
---|
| 3037 | import os |
---|
| 3038 | from Scientific.IO.NetCDF import NetCDFFile |
---|
| 3039 | from Numeric import Float, zeros, sum, reshape, equal |
---|
| 3040 | |
---|
| 3041 | root = basename_in |
---|
| 3042 | inname = root + '.dem' |
---|
| 3043 | |
---|
| 3044 | #Open existing netcdf file to read |
---|
| 3045 | infile = NetCDFFile(inname, 'r') |
---|
| 3046 | if verbose: print 'Reading DEM from %s' %inname |
---|
| 3047 | |
---|
| 3048 | #Read metadata |
---|
| 3049 | ncols = infile.ncols[0] |
---|
| 3050 | nrows = infile.nrows[0] |
---|
| 3051 | xllcorner = infile.xllcorner[0] |
---|
| 3052 | yllcorner = infile.yllcorner[0] |
---|
| 3053 | cellsize = infile.cellsize[0] |
---|
| 3054 | NODATA_value = infile.NODATA_value[0] |
---|
| 3055 | zone = infile.zone[0] |
---|
| 3056 | false_easting = infile.false_easting[0] |
---|
| 3057 | false_northing = infile.false_northing[0] |
---|
| 3058 | projection = infile.projection |
---|
| 3059 | datum = infile.datum |
---|
| 3060 | units = infile.units |
---|
| 3061 | |
---|
| 3062 | dem_elevation = infile.variables['elevation'] |
---|
| 3063 | |
---|
| 3064 | #Get output file name |
---|
| 3065 | if basename_out == None: |
---|
| 3066 | outname = root + '_' + repr(cellsize_new) + '.dem' |
---|
| 3067 | else: |
---|
| 3068 | outname = basename_out + '.dem' |
---|
| 3069 | |
---|
| 3070 | if verbose: print 'Write decimated NetCDF file to %s' %outname |
---|
| 3071 | |
---|
| 3072 | #Determine some dimensions for decimated grid |
---|
| 3073 | (nrows_stencil, ncols_stencil) = stencil.shape |
---|
| 3074 | x_offset = ncols_stencil / 2 |
---|
| 3075 | y_offset = nrows_stencil / 2 |
---|
| 3076 | cellsize_ratio = int(cellsize_new / cellsize) |
---|
| 3077 | ncols_new = 1 + (ncols - ncols_stencil) / cellsize_ratio |
---|
| 3078 | nrows_new = 1 + (nrows - nrows_stencil) / cellsize_ratio |
---|
| 3079 | |
---|
| 3080 | #Open netcdf file for output |
---|
| 3081 | outfile = NetCDFFile(outname, 'w') |
---|
| 3082 | |
---|
| 3083 | #Create new file |
---|
| 3084 | outfile.institution = 'Geoscience Australia' |
---|
| 3085 | outfile.description = 'NetCDF DEM format for compact and portable storage ' +\ |
---|
| 3086 | 'of spatial point data' |
---|
| 3087 | #Georeferencing |
---|
| 3088 | outfile.zone = zone |
---|
| 3089 | outfile.projection = projection |
---|
| 3090 | outfile.datum = datum |
---|
| 3091 | outfile.units = units |
---|
| 3092 | |
---|
| 3093 | outfile.cellsize = cellsize_new |
---|
| 3094 | outfile.NODATA_value = NODATA_value |
---|
| 3095 | outfile.false_easting = false_easting |
---|
| 3096 | outfile.false_northing = false_northing |
---|
| 3097 | |
---|
| 3098 | outfile.xllcorner = xllcorner + (x_offset * cellsize) |
---|
| 3099 | outfile.yllcorner = yllcorner + (y_offset * cellsize) |
---|
| 3100 | outfile.ncols = ncols_new |
---|
| 3101 | outfile.nrows = nrows_new |
---|
| 3102 | |
---|
| 3103 | # dimension definition |
---|
| 3104 | outfile.createDimension('number_of_points', nrows_new*ncols_new) |
---|
| 3105 | |
---|
| 3106 | # variable definition |
---|
| 3107 | outfile.createVariable('elevation', Float, ('number_of_points',)) |
---|
| 3108 | |
---|
| 3109 | # Get handle to the variable |
---|
| 3110 | elevation = outfile.variables['elevation'] |
---|
| 3111 | |
---|
| 3112 | dem_elevation_r = reshape(dem_elevation, (nrows, ncols)) |
---|
| 3113 | |
---|
| 3114 | #Store data |
---|
| 3115 | global_index = 0 |
---|
| 3116 | for i in range(nrows_new): |
---|
| 3117 | if verbose: print 'Processing row %d of %d' %(i, nrows_new) |
---|
| 3118 | lower_index = global_index |
---|
| 3119 | telev = zeros(ncols_new, Float) |
---|
| 3120 | local_index = 0 |
---|
| 3121 | trow = i * cellsize_ratio |
---|
| 3122 | |
---|
| 3123 | for j in range(ncols_new): |
---|
| 3124 | tcol = j * cellsize_ratio |
---|
| 3125 | tmp = dem_elevation_r[trow:trow+nrows_stencil, tcol:tcol+ncols_stencil] |
---|
| 3126 | |
---|
| 3127 | #if dem contains 1 or more NODATA_values set value in |
---|
| 3128 | #decimated dem to NODATA_value, else compute decimated |
---|
| 3129 | #value using stencil |
---|
| 3130 | if sum(sum(equal(tmp, NODATA_value))) > 0: |
---|
| 3131 | telev[local_index] = NODATA_value |
---|
| 3132 | else: |
---|
| 3133 | telev[local_index] = sum(sum(tmp * stencil)) |
---|
| 3134 | |
---|
| 3135 | global_index += 1 |
---|
| 3136 | local_index += 1 |
---|
| 3137 | |
---|
| 3138 | upper_index = global_index |
---|
| 3139 | |
---|
| 3140 | elevation[lower_index:upper_index] = telev |
---|
| 3141 | |
---|
| 3142 | assert global_index == nrows_new*ncols_new, 'index not equal to number of points' |
---|
| 3143 | |
---|
| 3144 | infile.close() |
---|
| 3145 | outfile.close() |
---|
| 3146 | |
---|
| 3147 | |
---|
| 3148 | |
---|
| 3149 | def sww2asc_obsolete(basename_in, basename_out = None, |
---|
| 3150 | quantity = None, |
---|
| 3151 | timestep = None, |
---|
| 3152 | reduction = None, |
---|
| 3153 | cellsize = 10, |
---|
| 3154 | verbose = False, |
---|
| 3155 | origin = None): |
---|
| 3156 | """Read SWW file and convert to Digitial Elevation model format (.asc) |
---|
| 3157 | |
---|
| 3158 | Example: |
---|
| 3159 | |
---|
| 3160 | ncols 3121 |
---|
| 3161 | nrows 1800 |
---|
| 3162 | xllcorner 722000 |
---|
| 3163 | yllcorner 5893000 |
---|
| 3164 | cellsize 25 |
---|
| 3165 | NODATA_value -9999 |
---|
| 3166 | 138.3698 137.4194 136.5062 135.5558 .......... |
---|
| 3167 | |
---|
| 3168 | Also write accompanying file with same basename_in but extension .prj |
---|
| 3169 | used to fix the UTM zone, datum, false northings and eastings. |
---|
| 3170 | |
---|
| 3171 | The prj format is assumed to be as |
---|
| 3172 | |
---|
| 3173 | Projection UTM |
---|
| 3174 | Zone 56 |
---|
| 3175 | Datum WGS84 |
---|
| 3176 | Zunits NO |
---|
| 3177 | Units METERS |
---|
| 3178 | Spheroid WGS84 |
---|
| 3179 | Xshift 0.0000000000 |
---|
| 3180 | Yshift 10000000.0000000000 |
---|
| 3181 | Parameters |
---|
| 3182 | |
---|
| 3183 | |
---|
| 3184 | if quantity is given, out values from quantity otherwise default to |
---|
| 3185 | elevation |
---|
| 3186 | |
---|
| 3187 | if timestep (an index) is given, output quantity at that timestep |
---|
| 3188 | |
---|
| 3189 | if reduction is given use that to reduce quantity over all timesteps. |
---|
| 3190 | |
---|
| 3191 | """ |
---|
| 3192 | from Numeric import array, Float, concatenate, NewAxis, zeros,\ |
---|
| 3193 | sometrue |
---|
| 3194 | from utilities.polygon import inside_polygon |
---|
| 3195 | |
---|
| 3196 | #FIXME: Should be variable |
---|
| 3197 | datum = 'WGS84' |
---|
| 3198 | false_easting = 500000 |
---|
| 3199 | false_northing = 10000000 |
---|
| 3200 | |
---|
| 3201 | if quantity is None: |
---|
| 3202 | quantity = 'elevation' |
---|
| 3203 | |
---|
| 3204 | if reduction is None: |
---|
| 3205 | reduction = max |
---|
| 3206 | |
---|
| 3207 | if basename_out is None: |
---|
| 3208 | basename_out = basename_in + '_%s' %quantity |
---|
| 3209 | |
---|
| 3210 | swwfile = basename_in + '.sww' |
---|
| 3211 | ascfile = basename_out + '.asc' |
---|
| 3212 | prjfile = basename_out + '.prj' |
---|
| 3213 | |
---|
| 3214 | |
---|
| 3215 | if verbose: print 'Reading from %s' %swwfile |
---|
| 3216 | #Read sww file |
---|
| 3217 | from Scientific.IO.NetCDF import NetCDFFile |
---|
| 3218 | fid = NetCDFFile(swwfile) |
---|
| 3219 | |
---|
| 3220 | #Get extent and reference |
---|
| 3221 | x = fid.variables['x'][:] |
---|
| 3222 | y = fid.variables['y'][:] |
---|
| 3223 | volumes = fid.variables['volumes'][:] |
---|
| 3224 | |
---|
| 3225 | ymin = min(y); ymax = max(y) |
---|
| 3226 | xmin = min(x); xmax = max(x) |
---|
| 3227 | |
---|
| 3228 | number_of_timesteps = fid.dimensions['number_of_timesteps'] |
---|
| 3229 | number_of_points = fid.dimensions['number_of_points'] |
---|
| 3230 | if origin is None: |
---|
| 3231 | |
---|
| 3232 | #Get geo_reference |
---|
| 3233 | #sww files don't have to have a geo_ref |
---|
| 3234 | try: |
---|
| 3235 | geo_reference = Geo_reference(NetCDFObject=fid) |
---|
| 3236 | except AttributeError, e: |
---|
| 3237 | geo_reference = Geo_reference() #Default georef object |
---|
| 3238 | |
---|
| 3239 | xllcorner = geo_reference.get_xllcorner() |
---|
| 3240 | yllcorner = geo_reference.get_yllcorner() |
---|
| 3241 | zone = geo_reference.get_zone() |
---|
| 3242 | else: |
---|
| 3243 | zone = origin[0] |
---|
| 3244 | xllcorner = origin[1] |
---|
| 3245 | yllcorner = origin[2] |
---|
| 3246 | |
---|
| 3247 | |
---|
| 3248 | #Get quantity and reduce if applicable |
---|
| 3249 | if verbose: print 'Reading quantity %s' %quantity |
---|
| 3250 | |
---|
| 3251 | if quantity.lower() == 'depth': |
---|
| 3252 | q = fid.variables['stage'][:] - fid.variables['elevation'][:] |
---|
| 3253 | else: |
---|
| 3254 | q = fid.variables[quantity][:] |
---|
| 3255 | |
---|
| 3256 | |
---|
| 3257 | if len(q.shape) == 2: |
---|
| 3258 | if verbose: print 'Reducing quantity %s' %quantity |
---|
| 3259 | q_reduced = zeros( number_of_points, Float ) |
---|
| 3260 | |
---|
| 3261 | for k in range(number_of_points): |
---|
| 3262 | q_reduced[k] = reduction( q[:,k] ) |
---|
| 3263 | |
---|
| 3264 | q = q_reduced |
---|
| 3265 | |
---|
| 3266 | #Now q has dimension: number_of_points |
---|
| 3267 | |
---|
| 3268 | #Create grid and update xll/yll corner and x,y |
---|
| 3269 | if verbose: print 'Creating grid' |
---|
| 3270 | ncols = int((xmax-xmin)/cellsize)+1 |
---|
| 3271 | nrows = int((ymax-ymin)/cellsize)+1 |
---|
| 3272 | |
---|
| 3273 | newxllcorner = xmin+xllcorner |
---|
| 3274 | newyllcorner = ymin+yllcorner |
---|
| 3275 | |
---|
| 3276 | x = x+xllcorner-newxllcorner |
---|
| 3277 | y = y+yllcorner-newyllcorner |
---|
| 3278 | |
---|
| 3279 | vertex_points = concatenate ((x[:, NewAxis] ,y[:, NewAxis]), axis = 1) |
---|
| 3280 | assert len(vertex_points.shape) == 2 |
---|
| 3281 | |
---|
| 3282 | |
---|
| 3283 | from Numeric import zeros, Float |
---|
| 3284 | grid_points = zeros ( (ncols*nrows, 2), Float ) |
---|
| 3285 | |
---|
| 3286 | |
---|
| 3287 | for i in xrange(nrows): |
---|
| 3288 | yg = i*cellsize |
---|
| 3289 | for j in xrange(ncols): |
---|
| 3290 | xg = j*cellsize |
---|
| 3291 | k = i*ncols + j |
---|
| 3292 | |
---|
| 3293 | grid_points[k,0] = xg |
---|
| 3294 | grid_points[k,1] = yg |
---|
| 3295 | |
---|
| 3296 | #Interpolate |
---|
| 3297 | from least_squares import Interpolation |
---|
| 3298 | |
---|
| 3299 | |
---|
| 3300 | #FIXME: This should be done with precrop = True, otherwise it'll |
---|
| 3301 | #take forever. With expand_search set to False, some grid points might |
---|
| 3302 | #miss out.... |
---|
| 3303 | |
---|
| 3304 | interp = Interpolation(vertex_points, volumes, grid_points, alpha=0.0, |
---|
| 3305 | precrop = False, expand_search = True, |
---|
| 3306 | verbose = verbose) |
---|
| 3307 | |
---|
| 3308 | #Interpolate using quantity values |
---|
| 3309 | if verbose: print 'Interpolating' |
---|
| 3310 | grid_values = interp.interpolate(q).flat |
---|
| 3311 | |
---|
| 3312 | #Write |
---|
| 3313 | #Write prj file |
---|
| 3314 | if verbose: print 'Writing %s' %prjfile |
---|
| 3315 | prjid = open(prjfile, 'w') |
---|
| 3316 | prjid.write('Projection %s\n' %'UTM') |
---|
| 3317 | prjid.write('Zone %d\n' %zone) |
---|
| 3318 | prjid.write('Datum %s\n' %datum) |
---|
| 3319 | prjid.write('Zunits NO\n') |
---|
| 3320 | prjid.write('Units METERS\n') |
---|
| 3321 | prjid.write('Spheroid %s\n' %datum) |
---|
| 3322 | prjid.write('Xshift %d\n' %false_easting) |
---|
| 3323 | prjid.write('Yshift %d\n' %false_northing) |
---|
| 3324 | prjid.write('Parameters\n') |
---|
| 3325 | prjid.close() |
---|
| 3326 | |
---|
| 3327 | |
---|
| 3328 | |
---|
| 3329 | if verbose: print 'Writing %s' %ascfile |
---|
| 3330 | NODATA_value = -9999 |
---|
| 3331 | |
---|
| 3332 | ascid = open(ascfile, 'w') |
---|
| 3333 | |
---|
| 3334 | ascid.write('ncols %d\n' %ncols) |
---|
| 3335 | ascid.write('nrows %d\n' %nrows) |
---|
| 3336 | ascid.write('xllcorner %d\n' %newxllcorner) |
---|
| 3337 | ascid.write('yllcorner %d\n' %newyllcorner) |
---|
| 3338 | ascid.write('cellsize %f\n' %cellsize) |
---|
| 3339 | ascid.write('NODATA_value %d\n' %NODATA_value) |
---|
| 3340 | |
---|
| 3341 | |
---|
| 3342 | #Get bounding polygon from mesh |
---|
| 3343 | P = interp.mesh.get_boundary_polygon() |
---|
| 3344 | inside_indices = inside_polygon(grid_points, P) |
---|
| 3345 | |
---|
| 3346 | for i in range(nrows): |
---|
| 3347 | if verbose and i%((nrows+10)/10)==0: |
---|
| 3348 | print 'Doing row %d of %d' %(i, nrows) |
---|
| 3349 | |
---|
| 3350 | for j in range(ncols): |
---|
| 3351 | index = (nrows-i-1)*ncols+j |
---|
| 3352 | |
---|
| 3353 | if sometrue(inside_indices == index): |
---|
| 3354 | ascid.write('%f ' %grid_values[index]) |
---|
| 3355 | else: |
---|
| 3356 | ascid.write('%d ' %NODATA_value) |
---|
| 3357 | |
---|
| 3358 | ascid.write('\n') |
---|
| 3359 | |
---|
| 3360 | #Close |
---|
| 3361 | ascid.close() |
---|
| 3362 | fid.close() |
---|
| 3363 | |
---|
| 3364 | #******************** |
---|
| 3365 | #*** END OF OBSOLETE FUNCTIONS |
---|
| 3366 | #*************** |
---|
| 3367 | |
---|
| 3368 | |
---|
| 3369 | def tsh2sww(filename, verbose=False): #test_tsh2sww |
---|
| 3370 | """ |
---|
| 3371 | to check if a tsh/msh file 'looks' good. |
---|
| 3372 | """ |
---|
| 3373 | |
---|
| 3374 | #FIXME Move to data_manager |
---|
| 3375 | from shallow_water import Domain |
---|
| 3376 | from pmesh2domain import pmesh_to_domain_instance |
---|
| 3377 | import time, os |
---|
| 3378 | from data_manager import get_dataobject |
---|
| 3379 | from os import sep, path |
---|
| 3380 | from util import mean |
---|
| 3381 | |
---|
| 3382 | if verbose == True:print 'Creating domain from', filename |
---|
| 3383 | domain = pmesh_to_domain_instance(filename, Domain) |
---|
| 3384 | if verbose == True:print "Number of triangles = ", len(domain) |
---|
| 3385 | |
---|
| 3386 | domain.smooth = True |
---|
| 3387 | domain.format = 'sww' #Native netcdf visualisation format |
---|
| 3388 | file_path, filename = path.split(filename) |
---|
| 3389 | filename, ext = path.splitext(filename) |
---|
| 3390 | domain.filename = filename |
---|
| 3391 | domain.reduction = mean |
---|
| 3392 | if verbose == True:print "file_path",file_path |
---|
| 3393 | if file_path == "":file_path = "." |
---|
| 3394 | domain.set_datadir(file_path) |
---|
| 3395 | |
---|
| 3396 | if verbose == True: |
---|
| 3397 | print "Output written to " + domain.get_datadir() + sep + \ |
---|
| 3398 | domain.filename + "." + domain.format |
---|
| 3399 | sww = get_dataobject(domain) |
---|
| 3400 | sww.store_connectivity() |
---|
| 3401 | sww.store_timestep('stage') |
---|
| 3402 | |
---|
| 3403 | |
---|
| 3404 | def asc_csiro2sww(bath_dir, |
---|
| 3405 | elevation_dir, |
---|
| 3406 | ucur_dir, |
---|
| 3407 | vcur_dir, |
---|
| 3408 | sww_file, |
---|
| 3409 | minlat = None, maxlat = None, |
---|
| 3410 | minlon = None, maxlon = None, |
---|
| 3411 | zscale=1, |
---|
| 3412 | mean_stage = 0, |
---|
| 3413 | fail_on_NaN = True, |
---|
| 3414 | elevation_NaN_filler = 0, |
---|
| 3415 | bath_prefix='ba', |
---|
| 3416 | elevation_prefix='el', |
---|
| 3417 | verbose=False): |
---|
| 3418 | """ |
---|
| 3419 | Produce an sww boundary file, from esri ascii data from CSIRO. |
---|
| 3420 | |
---|
| 3421 | Also convert latitude and longitude to UTM. All coordinates are |
---|
| 3422 | assumed to be given in the GDA94 datum. |
---|
| 3423 | |
---|
| 3424 | assume: |
---|
| 3425 | All files are in esri ascii format |
---|
| 3426 | |
---|
| 3427 | 4 types of information |
---|
| 3428 | bathymetry |
---|
| 3429 | elevation |
---|
| 3430 | u velocity |
---|
| 3431 | v velocity |
---|
| 3432 | |
---|
| 3433 | Assumptions |
---|
| 3434 | The metadata of all the files is the same |
---|
| 3435 | Each type is in a seperate directory |
---|
| 3436 | One bath file with extention .000 |
---|
| 3437 | The time period is less than 24hrs and uniform. |
---|
| 3438 | """ |
---|
| 3439 | from Scientific.IO.NetCDF import NetCDFFile |
---|
| 3440 | |
---|
| 3441 | from coordinate_transforms.redfearn import redfearn |
---|
| 3442 | |
---|
| 3443 | precision = Float # So if we want to change the precision its done here |
---|
| 3444 | |
---|
| 3445 | # go in to the bath dir and load the only file, |
---|
| 3446 | bath_files = os.listdir(bath_dir) |
---|
| 3447 | #print "bath_files",bath_files |
---|
| 3448 | |
---|
| 3449 | #fixme: make more general? |
---|
| 3450 | bath_file = bath_files[0] |
---|
| 3451 | bath_dir_file = bath_dir + os.sep + bath_file |
---|
| 3452 | bath_metadata,bath_grid = _read_asc(bath_dir_file) |
---|
| 3453 | #print "bath_metadata",bath_metadata |
---|
| 3454 | #print "bath_grid",bath_grid |
---|
| 3455 | |
---|
| 3456 | #Use the date.time of the bath file as a basis for |
---|
| 3457 | #the start time for other files |
---|
| 3458 | base_start = bath_file[-12:] |
---|
| 3459 | |
---|
| 3460 | #go into the elevation dir and load the 000 file |
---|
| 3461 | elevation_dir_file = elevation_dir + os.sep + elevation_prefix \ |
---|
| 3462 | + base_start |
---|
| 3463 | #elevation_metadata,elevation_grid = _read_asc(elevation_dir_file) |
---|
| 3464 | #print "elevation_dir_file",elevation_dir_file |
---|
| 3465 | #print "elevation_grid", elevation_grid |
---|
| 3466 | |
---|
| 3467 | elevation_files = os.listdir(elevation_dir) |
---|
| 3468 | ucur_files = os.listdir(ucur_dir) |
---|
| 3469 | vcur_files = os.listdir(vcur_dir) |
---|
| 3470 | |
---|
| 3471 | # the first elevation file should be the |
---|
| 3472 | # file with the same base name as the bath data |
---|
| 3473 | #print "elevation_files[0]",elevation_files[0] |
---|
| 3474 | #print "'el' + base_start",'el' + base_start |
---|
| 3475 | assert elevation_files[0] == 'el' + base_start |
---|
| 3476 | |
---|
| 3477 | #assert bath_metadata == elevation_metadata |
---|
| 3478 | |
---|
| 3479 | |
---|
| 3480 | |
---|
| 3481 | number_of_latitudes = bath_grid.shape[0] |
---|
| 3482 | number_of_longitudes = bath_grid.shape[1] |
---|
| 3483 | #number_of_times = len(os.listdir(elevation_dir)) |
---|
| 3484 | #number_of_points = number_of_latitudes*number_of_longitudes |
---|
| 3485 | number_of_volumes = (number_of_latitudes-1)*(number_of_longitudes-1)*2 |
---|
| 3486 | |
---|
| 3487 | longitudes = [bath_metadata['xllcorner']+x*bath_metadata['cellsize'] \ |
---|
| 3488 | for x in range(number_of_longitudes)] |
---|
| 3489 | latitudes = [bath_metadata['yllcorner']+y*bath_metadata['cellsize'] \ |
---|
| 3490 | for y in range(number_of_latitudes)] |
---|
| 3491 | |
---|
| 3492 | # reverse order of lat, so the fist lat represents the first grid row |
---|
| 3493 | latitudes.reverse() |
---|
| 3494 | |
---|
| 3495 | #print "latitudes - before _get_min_max_indexes",latitudes |
---|
| 3496 | kmin, kmax, lmin, lmax = _get_min_max_indexes(latitudes,longitudes, |
---|
| 3497 | minlat=minlat, maxlat=maxlat, |
---|
| 3498 | minlon=minlon, maxlon=maxlon) |
---|
| 3499 | |
---|
| 3500 | |
---|
| 3501 | bath_grid = bath_grid[kmin:kmax,lmin:lmax] |
---|
| 3502 | #print "bath_grid",bath_grid |
---|
| 3503 | latitudes = latitudes[kmin:kmax] |
---|
| 3504 | longitudes = longitudes[lmin:lmax] |
---|
| 3505 | number_of_latitudes = len(latitudes) |
---|
| 3506 | number_of_longitudes = len(longitudes) |
---|
| 3507 | number_of_times = len(os.listdir(elevation_dir)) |
---|
| 3508 | number_of_points = number_of_latitudes*number_of_longitudes |
---|
| 3509 | number_of_volumes = (number_of_latitudes-1)*(number_of_longitudes-1)*2 |
---|
| 3510 | #print "number_of_points",number_of_points |
---|
| 3511 | |
---|
| 3512 | #Work out the times |
---|
| 3513 | if len(elevation_files) > 1: |
---|
| 3514 | # Assume: The time period is less than 24hrs. |
---|
| 3515 | time_period = (int(elevation_files[1][-3:]) - \ |
---|
| 3516 | int(elevation_files[0][-3:]))*60*60 |
---|
| 3517 | times = [x*time_period for x in range(len(elevation_files))] |
---|
| 3518 | else: |
---|
| 3519 | times = [0.0] |
---|
| 3520 | #print "times", times |
---|
| 3521 | #print "number_of_latitudes", number_of_latitudes |
---|
| 3522 | #print "number_of_longitudes", number_of_longitudes |
---|
| 3523 | #print "number_of_times", number_of_times |
---|
| 3524 | #print "latitudes", latitudes |
---|
| 3525 | #print "longitudes", longitudes |
---|
| 3526 | |
---|
| 3527 | |
---|
| 3528 | if verbose: |
---|
| 3529 | print '------------------------------------------------' |
---|
| 3530 | print 'Statistics:' |
---|
| 3531 | print ' Extent (lat/lon):' |
---|
| 3532 | print ' lat in [%f, %f], len(lat) == %d'\ |
---|
| 3533 | %(min(latitudes), max(latitudes), |
---|
| 3534 | len(latitudes)) |
---|
| 3535 | print ' lon in [%f, %f], len(lon) == %d'\ |
---|
| 3536 | %(min(longitudes), max(longitudes), |
---|
| 3537 | len(longitudes)) |
---|
| 3538 | print ' t in [%f, %f], len(t) == %d'\ |
---|
| 3539 | %(min(times), max(times), len(times)) |
---|
| 3540 | |
---|
| 3541 | ######### WRITE THE SWW FILE ############# |
---|
| 3542 | # NetCDF file definition |
---|
| 3543 | outfile = NetCDFFile(sww_file, 'w') |
---|
| 3544 | |
---|
| 3545 | #Create new file |
---|
| 3546 | outfile.institution = 'Geoscience Australia' |
---|
| 3547 | outfile.description = 'Converted from XXX' |
---|
| 3548 | |
---|
| 3549 | |
---|
| 3550 | #For sww compatibility |
---|
| 3551 | outfile.smoothing = 'Yes' |
---|
| 3552 | outfile.order = 1 |
---|
| 3553 | |
---|
| 3554 | #Start time in seconds since the epoch (midnight 1/1/1970) |
---|
| 3555 | outfile.starttime = starttime = times[0] |
---|
| 3556 | |
---|
| 3557 | |
---|
| 3558 | # dimension definitions |
---|
| 3559 | outfile.createDimension('number_of_volumes', number_of_volumes) |
---|
| 3560 | |
---|
| 3561 | outfile.createDimension('number_of_vertices', 3) |
---|
| 3562 | outfile.createDimension('number_of_points', number_of_points) |
---|
| 3563 | outfile.createDimension('number_of_timesteps', number_of_times) |
---|
| 3564 | |
---|
| 3565 | # variable definitions |
---|
| 3566 | outfile.createVariable('x', precision, ('number_of_points',)) |
---|
| 3567 | outfile.createVariable('y', precision, ('number_of_points',)) |
---|
| 3568 | outfile.createVariable('elevation', precision, ('number_of_points',)) |
---|
| 3569 | |
---|
| 3570 | #FIXME: Backwards compatibility |
---|
| 3571 | outfile.createVariable('z', precision, ('number_of_points',)) |
---|
| 3572 | ################################# |
---|
| 3573 | |
---|
| 3574 | outfile.createVariable('volumes', Int, ('number_of_volumes', |
---|
| 3575 | 'number_of_vertices')) |
---|
| 3576 | |
---|
| 3577 | outfile.createVariable('time', precision, |
---|
| 3578 | ('number_of_timesteps',)) |
---|
| 3579 | |
---|
| 3580 | outfile.createVariable('stage', precision, |
---|
| 3581 | ('number_of_timesteps', |
---|
| 3582 | 'number_of_points')) |
---|
| 3583 | |
---|
| 3584 | outfile.createVariable('xmomentum', precision, |
---|
| 3585 | ('number_of_timesteps', |
---|
| 3586 | 'number_of_points')) |
---|
| 3587 | |
---|
| 3588 | outfile.createVariable('ymomentum', precision, |
---|
| 3589 | ('number_of_timesteps', |
---|
| 3590 | 'number_of_points')) |
---|
| 3591 | |
---|
| 3592 | #Store |
---|
| 3593 | from coordinate_transforms.redfearn import redfearn |
---|
| 3594 | x = zeros(number_of_points, Float) #Easting |
---|
| 3595 | y = zeros(number_of_points, Float) #Northing |
---|
| 3596 | |
---|
| 3597 | if verbose: print 'Making triangular grid' |
---|
| 3598 | #Get zone of 1st point. |
---|
| 3599 | refzone, _, _ = redfearn(latitudes[0],longitudes[0]) |
---|
| 3600 | |
---|
| 3601 | vertices = {} |
---|
| 3602 | i = 0 |
---|
| 3603 | for k, lat in enumerate(latitudes): |
---|
| 3604 | for l, lon in enumerate(longitudes): |
---|
| 3605 | |
---|
| 3606 | vertices[l,k] = i |
---|
| 3607 | |
---|
| 3608 | zone, easting, northing = redfearn(lat,lon) |
---|
| 3609 | |
---|
| 3610 | msg = 'Zone boundary crossed at longitude =', lon |
---|
| 3611 | #assert zone == refzone, msg |
---|
| 3612 | #print '%7.2f %7.2f %8.2f %8.2f' %(lon, lat, easting, northing) |
---|
| 3613 | x[i] = easting |
---|
| 3614 | y[i] = northing |
---|
| 3615 | i += 1 |
---|
| 3616 | |
---|
| 3617 | |
---|
| 3618 | #Construct 2 triangles per 'rectangular' element |
---|
| 3619 | volumes = [] |
---|
| 3620 | for l in range(number_of_longitudes-1): #X direction |
---|
| 3621 | for k in range(number_of_latitudes-1): #Y direction |
---|
| 3622 | v1 = vertices[l,k+1] |
---|
| 3623 | v2 = vertices[l,k] |
---|
| 3624 | v3 = vertices[l+1,k+1] |
---|
| 3625 | v4 = vertices[l+1,k] |
---|
| 3626 | |
---|
| 3627 | #Note, this is different to the ferrit2sww code |
---|
| 3628 | #since the order of the lats is reversed. |
---|
| 3629 | volumes.append([v1,v3,v2]) #Upper element |
---|
| 3630 | volumes.append([v4,v2,v3]) #Lower element |
---|
| 3631 | |
---|
| 3632 | volumes = array(volumes) |
---|
| 3633 | |
---|
| 3634 | geo_ref = Geo_reference(refzone,min(x),min(y)) |
---|
| 3635 | geo_ref.write_NetCDF(outfile) |
---|
| 3636 | |
---|
| 3637 | # This will put the geo ref in the middle |
---|
| 3638 | #geo_ref = Geo_reference(refzone,(max(x)+min(x))/2.0,(max(x)+min(y))/2.) |
---|
| 3639 | |
---|
| 3640 | |
---|
| 3641 | if verbose: |
---|
| 3642 | print '------------------------------------------------' |
---|
| 3643 | print 'More Statistics:' |
---|
| 3644 | print ' Extent (/lon):' |
---|
| 3645 | print ' x in [%f, %f], len(lat) == %d'\ |
---|
| 3646 | %(min(x), max(x), |
---|
| 3647 | len(x)) |
---|
| 3648 | print ' y in [%f, %f], len(lon) == %d'\ |
---|
| 3649 | %(min(y), max(y), |
---|
| 3650 | len(y)) |
---|
| 3651 | print 'geo_ref: ',geo_ref |
---|
| 3652 | |
---|
| 3653 | z = resize(bath_grid,outfile.variables['z'][:].shape) |
---|
| 3654 | outfile.variables['x'][:] = x - geo_ref.get_xllcorner() |
---|
| 3655 | outfile.variables['y'][:] = y - geo_ref.get_yllcorner() |
---|
| 3656 | outfile.variables['z'][:] = z |
---|
| 3657 | outfile.variables['elevation'][:] = z #FIXME HACK |
---|
| 3658 | outfile.variables['volumes'][:] = volumes.astype(Int32) #On Opteron 64 |
---|
| 3659 | |
---|
| 3660 | # do this to create an ok sww file. |
---|
| 3661 | #outfile.variables['time'][:] = [0] #Store time relative |
---|
| 3662 | #outfile.variables['stage'] = z |
---|
| 3663 | # put the next line up in the code after outfile.order = 1 |
---|
| 3664 | #number_of_times = 1 |
---|
| 3665 | |
---|
| 3666 | stage = outfile.variables['stage'] |
---|
| 3667 | xmomentum = outfile.variables['xmomentum'] |
---|
| 3668 | ymomentum = outfile.variables['ymomentum'] |
---|
| 3669 | |
---|
| 3670 | outfile.variables['time'][:] = times #Store time relative |
---|
| 3671 | |
---|
| 3672 | if verbose: print 'Converting quantities' |
---|
| 3673 | n = number_of_times |
---|
| 3674 | for j in range(number_of_times): |
---|
| 3675 | # load in files |
---|
| 3676 | elevation_meta, elevation_grid = \ |
---|
| 3677 | _read_asc(elevation_dir + os.sep + elevation_files[j]) |
---|
| 3678 | |
---|
| 3679 | _, u_momentum_grid = _read_asc(ucur_dir + os.sep + ucur_files[j]) |
---|
| 3680 | _, v_momentum_grid = _read_asc(vcur_dir + os.sep + vcur_files[j]) |
---|
| 3681 | |
---|
| 3682 | #print "elevation_grid",elevation_grid |
---|
| 3683 | #cut matrix to desired size |
---|
| 3684 | elevation_grid = elevation_grid[kmin:kmax,lmin:lmax] |
---|
| 3685 | u_momentum_grid = u_momentum_grid[kmin:kmax,lmin:lmax] |
---|
| 3686 | v_momentum_grid = v_momentum_grid[kmin:kmax,lmin:lmax] |
---|
| 3687 | #print "elevation_grid",elevation_grid |
---|
| 3688 | # handle missing values |
---|
| 3689 | missing = (elevation_grid == elevation_meta['NODATA_value']) |
---|
| 3690 | if sometrue (missing): |
---|
| 3691 | if fail_on_NaN: |
---|
| 3692 | msg = 'File %s contains missing values'\ |
---|
| 3693 | %(elevation_files[j]) |
---|
| 3694 | raise msg |
---|
| 3695 | else: |
---|
| 3696 | elevation_grid = elevation_grid*(missing==0) + \ |
---|
| 3697 | missing*elevation_NaN_filler |
---|
| 3698 | |
---|
| 3699 | |
---|
| 3700 | if verbose and j%((n+10)/10)==0: print ' Doing %d of %d' %(j, n) |
---|
| 3701 | i = 0 |
---|
| 3702 | for k in range(number_of_latitudes): #Y direction |
---|
| 3703 | for l in range(number_of_longitudes): #X direction |
---|
| 3704 | w = zscale*elevation_grid[k,l] + mean_stage |
---|
| 3705 | stage[j,i] = w |
---|
| 3706 | h = w - z[i] |
---|
| 3707 | xmomentum[j,i] = u_momentum_grid[k,l]*h |
---|
| 3708 | ymomentum[j,i] = v_momentum_grid[k,l]*h |
---|
| 3709 | i += 1 |
---|
| 3710 | outfile.close() |
---|
| 3711 | |
---|
| 3712 | def _get_min_max_indexes(latitudes,longitudes, |
---|
| 3713 | minlat=None, maxlat=None, |
---|
| 3714 | minlon=None, maxlon=None): |
---|
| 3715 | """ |
---|
| 3716 | return max, min indexes (for slicing) of the lat and long arrays to cover the area |
---|
| 3717 | specified with min/max lat/long |
---|
| 3718 | |
---|
| 3719 | Think of the latitudes and longitudes describing a 2d surface. |
---|
| 3720 | The area returned is, if possible, just big enough to cover the |
---|
| 3721 | inputed max/min area. (This will not be possible if the max/min area |
---|
| 3722 | has a section outside of the latitudes/longitudes area.) |
---|
| 3723 | |
---|
| 3724 | assume latitudess & longitudes are sorted, |
---|
| 3725 | long - from low to high (west to east, eg 148 - 151) |
---|
| 3726 | lat - from high to low (north to south, eg -35 - -38) |
---|
| 3727 | """ |
---|
| 3728 | |
---|
| 3729 | # reverse order of lat, so it's in ascending order |
---|
| 3730 | latitudes.reverse() |
---|
| 3731 | largest_lat_index = len(latitudes)-1 |
---|
| 3732 | #Cut out a smaller extent. |
---|
| 3733 | if minlat == None: |
---|
| 3734 | lat_min_index = 0 |
---|
| 3735 | else: |
---|
| 3736 | lat_min_index = searchsorted(latitudes, minlat)-1 |
---|
| 3737 | if lat_min_index <0: |
---|
| 3738 | lat_min_index = 0 |
---|
| 3739 | |
---|
| 3740 | |
---|
| 3741 | if maxlat == None: |
---|
| 3742 | lat_max_index = largest_lat_index #len(latitudes) |
---|
| 3743 | else: |
---|
| 3744 | lat_max_index = searchsorted(latitudes, maxlat) |
---|
| 3745 | if lat_max_index > largest_lat_index: |
---|
| 3746 | lat_max_index = largest_lat_index |
---|
| 3747 | |
---|
| 3748 | if minlon == None: |
---|
| 3749 | lon_min_index = 0 |
---|
| 3750 | else: |
---|
| 3751 | lon_min_index = searchsorted(longitudes, minlon)-1 |
---|
| 3752 | if lon_min_index <0: |
---|
| 3753 | lon_min_index = 0 |
---|
| 3754 | |
---|
| 3755 | if maxlon == None: |
---|
| 3756 | lon_max_index = len(longitudes) |
---|
| 3757 | else: |
---|
| 3758 | lon_max_index = searchsorted(longitudes, maxlon) |
---|
| 3759 | |
---|
| 3760 | #Take into account that the latitude list was reversed |
---|
| 3761 | latitudes.reverse() # Python passes by reference, need to swap back |
---|
| 3762 | lat_min_index, lat_max_index = largest_lat_index - lat_max_index , \ |
---|
| 3763 | largest_lat_index - lat_min_index |
---|
| 3764 | lat_max_index = lat_max_index + 1 # taking into account how slicing works |
---|
| 3765 | lon_max_index = lon_max_index + 1 # taking into account how slicing works |
---|
| 3766 | |
---|
| 3767 | return lat_min_index, lat_max_index, lon_min_index, lon_max_index |
---|
| 3768 | |
---|
| 3769 | |
---|
| 3770 | def _read_asc(filename, verbose=False): |
---|
| 3771 | """Read esri file from the following ASCII format (.asc) |
---|
| 3772 | |
---|
| 3773 | Example: |
---|
| 3774 | |
---|
| 3775 | ncols 3121 |
---|
| 3776 | nrows 1800 |
---|
| 3777 | xllcorner 722000 |
---|
| 3778 | yllcorner 5893000 |
---|
| 3779 | cellsize 25 |
---|
| 3780 | NODATA_value -9999 |
---|
| 3781 | 138.3698 137.4194 136.5062 135.5558 .......... |
---|
| 3782 | |
---|
| 3783 | """ |
---|
| 3784 | |
---|
| 3785 | datafile = open(filename) |
---|
| 3786 | |
---|
| 3787 | if verbose: print 'Reading DEM from %s' %(filename) |
---|
| 3788 | lines = datafile.readlines() |
---|
| 3789 | datafile.close() |
---|
| 3790 | |
---|
| 3791 | if verbose: print 'Got', len(lines), ' lines' |
---|
| 3792 | |
---|
| 3793 | ncols = int(lines.pop(0).split()[1].strip()) |
---|
| 3794 | nrows = int(lines.pop(0).split()[1].strip()) |
---|
| 3795 | xllcorner = float(lines.pop(0).split()[1].strip()) |
---|
| 3796 | yllcorner = float(lines.pop(0).split()[1].strip()) |
---|
| 3797 | cellsize = float(lines.pop(0).split()[1].strip()) |
---|
| 3798 | NODATA_value = float(lines.pop(0).split()[1].strip()) |
---|
| 3799 | |
---|
| 3800 | assert len(lines) == nrows |
---|
| 3801 | |
---|
| 3802 | #Store data |
---|
| 3803 | grid = [] |
---|
| 3804 | |
---|
| 3805 | n = len(lines) |
---|
| 3806 | for i, line in enumerate(lines): |
---|
| 3807 | cells = line.split() |
---|
| 3808 | assert len(cells) == ncols |
---|
| 3809 | grid.append(array([float(x) for x in cells])) |
---|
| 3810 | grid = array(grid) |
---|
| 3811 | |
---|
| 3812 | return {'xllcorner':xllcorner, |
---|
| 3813 | 'yllcorner':yllcorner, |
---|
| 3814 | 'cellsize':cellsize, |
---|
| 3815 | 'NODATA_value':NODATA_value}, grid |
---|