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