Changeset 4455
- Timestamp:
- May 16, 2007, 3:13:52 PM (18 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
anuga_core/source/anuga/shallow_water/data_manager.py
r4418 r4455 55 55 class DataTimeError(exceptions.Exception): pass 56 56 class DataDomainError(exceptions.Exception): pass 57 class NewQuantity(exceptions.Exception): pass 57 58 58 59 … … 66 67 67 68 from Numeric import concatenate, array, Float, Int, Int32, resize, sometrue, \ 68 searchsorted, zeros, allclose, around, reshape, transpose, sort 69 searchsorted, zeros, allclose, around, reshape, transpose, sort, \ 70 NewAxis, ArrayType 69 71 from Scientific.IO.NetCDF import NetCDFFile 70 72 … … 73 75 convert_from_latlon_to_utm 74 76 from anuga.coordinate_transforms.geo_reference import Geo_reference, \ 75 write_NetCDF_georeference 77 write_NetCDF_georeference, ensure_geo_reference 76 78 from anuga.geospatial_data.geospatial_data import Geospatial_data,\ 77 79 ensure_absolute … … 83 85 from anuga.abstract_2d_finite_volumes.pmesh2domain import \ 84 86 pmesh_to_domain_instance 87 from anuga.abstract_2d_finite_volumes.util import get_revision_number 85 88 86 89 def make_filename(s): … … 280 283 281 284 if mode == 'w': 282 283 #Create new file 284 fid.institution = 'Geoscience Australia' 285 fid.description = 'Output from anuga.abstract_2d_finite_volumes suitable for plotting' 286 287 if domain.smooth: 288 fid.smoothing = 'Yes' 289 else: 290 fid.smoothing = 'No' 291 292 fid.order = domain.default_order 293 285 description = 'Output from anuga.abstract_2d_finite_volumes suitable for plotting' 286 self.writer = Write_sww() 287 self.writer.header(fid, domain.starttime, 288 self.number_of_volumes, 289 self.domain.number_of_full_nodes, 290 description=description, 291 smoothing=domain.smooth, 292 order=domain.default_order) 294 293 if hasattr(domain, 'texture'): 295 fid.texture = domain.texture 296 #else: 297 # fid.texture = 'None' 298 299 #Reference point 300 #Start time in seconds since the epoch (midnight 1/1/1970) 301 #FIXME: Use Georef 302 fid.starttime = domain.starttime 303 304 # dimension definitions 305 fid.createDimension('number_of_volumes', self.number_of_volumes) 306 fid.createDimension('number_of_vertices', 3) 307 308 if domain.smooth is True: 309 #fid.createDimension('number_of_points', len(domain.vertexlist)) 310 fid.createDimension('number_of_points', self.domain.number_of_full_nodes) 311 312 # FIXME(Ole): This will cause sww files for paralle domains to 313 # have ghost nodes stored (but not used by triangles). 314 # To clean this up, we have to change get_vertex_values and friends in 315 # quantity.py (but I can't be bothered right now) 316 else: 317 fid.createDimension('number_of_points', 3*self.number_of_volumes) 318 319 fid.createDimension('number_of_timesteps', None) #extensible 320 321 # variable definitions 322 fid.createVariable('x', self.precision, ('number_of_points',)) 323 fid.createVariable('y', self.precision, ('number_of_points',)) 324 fid.createVariable('elevation', self.precision, ('number_of_points',)) 325 if domain.geo_reference is not None: 326 domain.geo_reference.write_NetCDF(fid) 327 328 #FIXME: Backwards compatibility 329 fid.createVariable('z', self.precision, ('number_of_points',)) 330 ################################# 331 332 fid.createVariable('volumes', Int, ('number_of_volumes', 333 'number_of_vertices')) 334 335 fid.createVariable('time', Float, # Always use full precision lest two timesteps 336 # close to each other may appear as the same step 337 ('number_of_timesteps',)) 338 339 fid.createVariable('stage', self.precision, 340 ('number_of_timesteps', 341 'number_of_points')) 342 343 fid.createVariable('xmomentum', self.precision, 344 ('number_of_timesteps', 345 'number_of_points')) 346 347 fid.createVariable('ymomentum', self.precision, 348 ('number_of_timesteps', 349 'number_of_points')) 294 fid.texture = domain.texture 295 # if domain.geo_reference is not None: 296 # domain.geo_reference.write_NetCDF(fid) 297 298 # fid.institution = 'Geoscience Australia' 299 # fid.description = 'Output from anuga.abstract_2d_finite_volumes suitable for plotting' 300 301 # if domain.smooth: 302 # fid.smoothing = 'Yes' 303 # else: 304 # fid.smoothing = 'No' 305 306 # fid.order = domain.default_order 307 308 # #Reference point 309 # #Start time in seconds since the epoch (midnight 1/1/1970) 310 # #FIXME: Use Georef 311 # fid.starttime = domain.starttime 312 313 # # dimension definitions 314 # fid.createDimension('number_of_volumes', self.number_of_volumes) 315 # fid.createDimension('number_of_vertices', 3) 316 317 # if domain.smooth is True: 318 # #fid.createDimension('number_of_points', len(domain.vertexlist)) 319 # fid.createDimension('number_of_points', self.domain.number_of_full_nodes) 320 321 # # FIXME(Ole): This will cause sww files for paralle domains to 322 # # have ghost nodes stored (but not used by triangles). 323 # # To clean this up, we have to change get_vertex_values and friends in 324 # # quantity.py (but I can't be bothered right now) 325 # else: 326 # fid.createDimension('number_of_points', 3*self.number_of_volumes) 327 328 # fid.createDimension('number_of_timesteps', None) #extensible 329 330 # # variable definitions 331 # fid.createVariable('x', self.precision, ('number_of_points',)) 332 # fid.createVariable('y', self.precision, ('number_of_points',)) 333 # fid.createVariable('elevation', self.precision, ('number_of_points',)) 334 # if domain.geo_reference is not None: 335 # domain.geo_reference.write_NetCDF(fid) 336 337 # #FIXME: Backwards compatibility 338 # fid.createVariable('z', self.precision, ('number_of_points',)) 339 # ################################# 340 341 # fid.createVariable('volumes', Int, ('number_of_volumes', 342 # 'number_of_vertices')) 343 344 # fid.createVariable('time', Float, # Always use full precision lest two timesteps 345 # # close to each other may appear as the same step 346 # ('number_of_timesteps',)) 347 348 # fid.createVariable('stage', self.precision, 349 # ('number_of_timesteps', 350 # 'number_of_points')) 351 352 # fid.createVariable('xmomentum', self.precision, 353 # ('number_of_timesteps', 354 # 'number_of_points')) 355 356 # fid.createVariable('ymomentum', self.precision, 357 # ('number_of_timesteps', 358 # 'number_of_points')) 350 359 351 360 #Close … … 381 390 precision=self.precision) 382 391 383 volumes[:] = V.astype(volumes.typecode()) 384 x[:] = X 385 y[:] = Y 386 z[:] = Z 387 388 #FIXME: Backwards compatibility 389 z = fid.variables['z'] 390 z[:] = Z 392 # 393 points = concatenate( (X[:,NewAxis],Y[:,NewAxis]), axis=1 ) 394 self.writer.triangulation(fid, points, 395 V.astype(volumes.typecode()), 396 Z, 397 points_georeference= \ 398 domain.geo_reference) 399 # 400 401 # volumes[:] = V.astype(volumes.typecode()) 402 # x[:] = X 403 # y[:] = Y 404 # z[:] = Z 405 406 # #FIXME: Backwards compatibility 407 # z = fid.variables['z'] 408 # z[:] = Z 391 409 ################################ 392 410 … … 405 423 406 424 from Numeric import choose 407 408 #Get NetCDF 425 426 #Get NetCDF 409 427 retries = 0 410 428 file_open = False … … 486 504 ymomentum = fid.variables['ymomentum'] 487 505 i = len(time) 488 489 #Store time490 time[i] = self.domain.time491 492 493 506 if type(names) not in [types.ListType, types.TupleType]: 494 507 names = [names] 495 508 496 for name in names: 497 # Get quantity 498 Q = domain.quantities[name] 499 A,V = Q.get_vertex_values(xy = False, 509 if 'stage' in names and 'xmomentum' in names and \ 510 'ymomentum' in names: 511 512 # Get stage 513 Q = domain.quantities['stage'] 514 A,_ = Q.get_vertex_values(xy = False, 515 precision = self.precision) 516 z = fid.variables['elevation'] 517 stage = choose(A-z[:] >= self.minimum_storable_height, 518 (z[:], A)) 519 520 # Get xmomentum 521 Q = domain.quantities['xmomentum'] 522 xmomentum, _ = Q.get_vertex_values(xy = False, 500 523 precision = self.precision) 501 502 503 #FIXME: Make this general (see below) 504 if name == 'stage': 505 z = fid.variables['elevation'] 506 #print z[:] 507 #print A-z[:] 508 509 510 A = choose(A-z[:] >= self.minimum_storable_height, 511 (z[:], A)) 512 stage[i,:] = A.astype(self.precision) 513 elif name == 'xmomentum': 514 xmomentum[i,:] = A.astype(self.precision) 515 elif name == 'ymomentum': 516 ymomentum[i,:] = A.astype(self.precision) 517 518 #As in.... 519 #eval( name + '[i,:] = A.astype(self.precision)' ) 520 #FIXME: But we need a UNIT test for that before refactoring 524 525 # Get ymomentum 526 Q = domain.quantities['ymomentum'] 527 ymomentum, _ = Q.get_vertex_values(xy = False, 528 precision = self.precision) 529 self.writer.quantities(fid, 530 time=self.domain.time, 531 precision=self.precision, 532 stage=stage, 533 xmomentum=xmomentum, 534 ymomentum=ymomentum) 535 else: 536 # This is producing a sww that is not standard. 537 #Store time 538 time[i] = self.domain.time 539 540 for name in names: 541 # Get quantity 542 Q = domain.quantities[name] 543 A,V = Q.get_vertex_values(xy = False, 544 precision = self.precision) 545 546 #FIXME: Make this general (see below) 547 if name == 'stage': 548 z = fid.variables['elevation'] 549 A = choose(A-z[:] >= self.minimum_storable_height, 550 (z[:], A)) 551 stage[i,:] = A.astype(self.precision) 552 elif name == 'xmomentum': 553 xmomentum[i,:] = A.astype(self.precision) 554 elif name == 'ymomentum': 555 ymomentum[i,:] = A.astype(self.precision) 556 557 #As in.... 558 #eval( name + '[i,:] = A.astype(self.precision)' ) 559 #FIXME: But we need a UNIT test for that before 560 # refactoring 521 561 522 562 … … 1471 1511 """ 1472 1512 1473 #FIXME: Can this be written feasibly using write_pts?1474 1475 1513 import os 1476 1514 from Scientific.IO.NetCDF import NetCDFFile … … 1554 1592 ptsname = basename_out + '.pts' 1555 1593 1556 #FIXME (DSG-ON): use loadASCII export_points_file 1557 if verbose: print 'Store to NetCDF file %s' %ptsname 1558 # NetCDF file definition 1559 outfile = NetCDFFile(ptsname, 'w') 1560 1561 #Create new file 1562 outfile.institution = 'Geoscience Australia' 1563 outfile.description = 'NetCDF pts format for compact and portable ' +\ 1564 'storage of spatial point data derived from HEC-RAS' 1565 1566 #Georeferencing 1567 outfile.zone = zone 1568 outfile.xllcorner = 0.0 1569 outfile.yllcorner = 0.0 1570 outfile.false_easting = 500000 #FIXME: Use defaults from e.g. config.py 1571 outfile.false_northing = 1000000 1572 1573 outfile.projection = projection 1574 outfile.datum = datum 1575 outfile.units = units 1576 1577 1578 outfile.createDimension('number_of_points', len(points)) 1579 outfile.createDimension('number_of_dimensions', 2) #This is 2d data 1580 1581 # variable definitions 1582 outfile.createVariable('points', Float, ('number_of_points', 1583 'number_of_dimensions')) 1584 outfile.createVariable('elevation', Float, ('number_of_points',)) 1585 1586 # Get handles to the variables 1587 outfile.variables['points'][:] = points 1588 outfile.variables['elevation'][:] = elevation 1589 1590 outfile.close() 1591 1594 geo_ref = Geo_reference(zone, 0, 0, datum, projection, units) 1595 geo = Geospatial_data(points, {"elevation":elevation}, 1596 verbose=verbose, geo_reference=geo_ref) 1597 geo.export_points_file(ptsname) 1592 1598 1593 1599 … … 2657 2663 #Create new file 2658 2664 starttime = times[0] 2659 write_sww_header(outfile, times, number_of_volumes, 2660 number_of_points, description=description) 2665 sww = Write_sww() 2666 sww.header(outfile, times, number_of_volumes, 2667 number_of_points, description=description, 2668 verbose=verbose) 2661 2669 2662 2670 #Store … … 2749 2757 x = outfile.variables['x'][:] 2750 2758 y = outfile.variables['y'][:] 2751 times = outfile.variables['time'][:]2752 2759 print '------------------------------------------------' 2753 2760 print 'Statistics of output file:' … … 3011 3018 other_quantities.remove('z') 3012 3019 other_quantities.remove('volumes') 3020 try: 3021 other_quantities.remove('stage_range') 3022 other_quantities.remove('xmomentum_range') 3023 other_quantities.remove('ymomentum_range') 3024 other_quantities.remove('elevation_range') 3025 except: 3026 pass 3027 3013 3028 3014 3029 conserved_quantities.remove('time') … … 3140 3155 3141 3156 def weed(coordinates,volumes,boundary = None): 3142 if type(coordinates)== 'array':3157 if type(coordinates)==ArrayType: 3143 3158 coordinates = coordinates.tolist() 3144 if type(volumes)== 'array':3159 if type(volumes)==ArrayType: 3145 3160 volumes = volumes.tolist() 3146 3161 … … 4367 4382 It will be the origin of the sww file. This shouldn't be used, 4368 4383 since all of anuga should be able to handle an arbitary origin. 4384 The mux point info is NOT relative to this origin. 4369 4385 4370 4386 … … 4437 4453 mesh = Mesh() 4438 4454 mesh.add_vertices(points_utm) 4439 mesh.auto_segment() 4455 mesh.auto_segment(smooth_indents=True, expand_pinch=True) 4456 # To try and avoid alpha shape 'hugging' too much 4457 mesh.auto_segment( mesh.shape.get_alpha()*1.1 ) 4440 4458 if hole_points_UTM is not None: 4441 4459 point = ensure_absolute(hole_points_UTM) … … 4485 4503 # For a different way of doing this, check out tsh2sww 4486 4504 # work out sww_times and the index range this covers 4487 write_sww_header(outfile, times, len(volumes), len(points_utm)) 4505 sww = Write_sww() 4506 sww.header(outfile, times, len(volumes), len(points_utm), 4507 verbose=verbose) 4488 4508 outfile.mean_stage = mean_stage 4489 4509 outfile.zscale = zscale 4490 write_sww_triangulation(outfile, points_utm, volumes, 4491 elevation, zone, origin=origin, verbose=verbose) 4510 4511 sww.triangulation(outfile, points_utm, volumes, 4512 elevation, zone, new_origin=origin, 4513 verbose=verbose) 4492 4514 4493 4515 if verbose: print 'Converting quantities' 4494 4516 j = 0 4495 4517 # Read in a time slice from each mux file and write it to the sww file 4496 for HA, UA, VAin map(None, mux['HA'], mux['UA'], mux['VA']):4518 for ha, ua, va in map(None, mux['HA'], mux['UA'], mux['VA']): 4497 4519 if j >= mux_times_start_i and j < mux_times_fin_i: 4498 write_sww_time_slice(outfile, HA, UA, VA, elevation, 4499 j - mux_times_start_i, 4500 mean_stage=mean_stage, zscale=zscale, 4501 verbose=verbose) 4520 stage = zscale*ha + mean_stage 4521 h = stage - elevation 4522 xmomentum = ua*h 4523 ymomentum = -1*va*h # -1 since in mux files south is positive. 4524 sww.quantities(outfile, 4525 slice_index=j - mux_times_start_i, 4526 verbose=verbose, 4527 stage=stage, 4528 xmomentum=xmomentum, 4529 ymomentum=ymomentum) 4502 4530 j += 1 4531 if verbose: sww.verbose_quantities(outfile) 4503 4532 outfile.close() 4504 4533 # 4505 4534 # Do some conversions while writing the sww file 4535 4536 4506 4537 def mux2sww_time(mux_times, mint, maxt): 4507 4538 """ … … 4522 4553 return mux_times_start_i, mux_times_fin_i 4523 4554 4524 def write_sww_header(outfile, times, number_of_volumes, number_of_points, description='Converted from XXX'): 4525 """ 4526 outfile - the name of the file that will be written 4527 times - A list of the time slice times 4528 number_of_volumes - the number of triangles 4529 """ 4530 times = ensure_numeric(times) 4531 4532 number_of_times = len(times) 4533 4534 outfile.institution = 'Geoscience Australia' 4535 outfile.description = 'Converted from XXX' 4536 4537 #For sww compatibility 4538 outfile.smoothing = 'Yes' 4539 outfile.order = 1 4540 4541 #Start time in seconds since the epoch (midnight 1/1/1970) 4542 if len(times) == 0: 4543 outfile.starttime = starttime = 0 4544 else: 4545 outfile.starttime = starttime = times[0] 4546 times = times - starttime #Store relative times 4547 4548 # dimension definitions 4549 outfile.createDimension('number_of_volumes', number_of_volumes) 4550 outfile.createDimension('number_of_vertices', 3) 4551 outfile.createDimension('number_of_points', number_of_points) 4552 outfile.createDimension('number_of_timesteps', number_of_times) 4553 4554 # variable definitions 4555 outfile.createVariable('x', precision, ('number_of_points',)) 4556 outfile.createVariable('y', precision, ('number_of_points',)) 4557 outfile.createVariable('elevation', precision, ('number_of_points',)) 4558 4559 #FIXME: Backwards compatibility 4560 outfile.createVariable('z', precision, ('number_of_points',)) 4561 ################################# 4562 4563 outfile.createVariable('volumes', Int, ('number_of_volumes', 4564 'number_of_vertices')) 4565 4566 outfile.createVariable('time', precision, 4567 ('number_of_timesteps',)) 4568 4569 outfile.createVariable('stage', precision, 4570 ('number_of_timesteps', 4571 'number_of_points')) 4572 4573 outfile.createVariable('xmomentum', precision, 4574 ('number_of_timesteps', 4575 'number_of_points')) 4576 4577 outfile.createVariable('ymomentum', precision, 4578 ('number_of_timesteps', 4579 'number_of_points')) 4580 4581 outfile.variables['time'][:] = times #Store time relative 4582 4583 4584 def write_sww_triangulation(outfile, points_utm, volumes, 4585 elevation, zone, origin=None, verbose=False): 4586 4587 number_of_points = len(points_utm) 4588 volumes = array(volumes) 4589 4590 if origin is None: 4591 origin = Geo_reference(zone,min(points_utm[:,0]),min(points_utm[:,1])) 4592 geo_ref = write_NetCDF_georeference(origin, outfile) 4593 4594 # if origin is not None: 4595 # if isinstance(origin, Geo_reference): 4596 # geo_ref = origin 4597 # else: 4598 # geo_ref = apply(Geo_reference,origin) 4599 # else: 4600 # geo_ref = Geo_reference(zone,min(points_utm[:,0]),min(points_utm[:,1])) 4601 # #geo_ref = Geo_reference(zone,0,0) 4602 # geo_ref.write_NetCDF(outfile) 4603 4604 4605 # This will put the geo ref in the middle 4606 #geo_ref = Geo_reference(refzone,(max(x)+min(x))/2.0,(max(x)+min(y))/2.) 4607 x = points_utm[:,0] 4608 y = points_utm[:,1] 4609 z = outfile.variables['z'][:] 4610 #FIXME the w.r.t time is lost now.. 4611 if verbose: 4555 4556 class Write_sww: 4557 from anuga.shallow_water.shallow_water_domain import Domain 4558 sww_quantities = Domain.conserved_quantities 4559 RANGE = '_range' 4560 def __init__(self): 4561 pass 4562 4563 def header(self,outfile, times, number_of_volumes, 4564 number_of_points, 4565 description='Converted from XXX', 4566 smoothing=True, 4567 order=1, verbose=False): 4568 """ 4569 outfile - the name of the file that will be written 4570 times - A list of the time slice times OR a start time 4571 Note, if a list is given the info will be made relative. 4572 number_of_volumes - the number of triangles 4573 """ 4574 4575 outfile.institution = 'Geoscience Australia' 4576 outfile.description = description 4577 4578 #For sww compatibility 4579 if smoothing is True: 4580 # Smoothing to be depreciated 4581 outfile.smoothing = 'Yes' 4582 outfile.vertices_are_stored_uniquely = 'False' 4583 else: 4584 # Smoothing to be depreciated 4585 outfile.smoothing = 'No' 4586 outfile.vertices_are_stored_uniquely = 'True' 4587 outfile.order = order 4588 4589 try: 4590 revision_number = get_revision_number() 4591 except: 4592 revision_number = None 4593 # writing a string so None can be written 4594 outfile.revision_number = str(revision_number) 4595 #times - A list or array of the time slice times OR a start time 4596 #times = ensure_numeric(times) 4597 #Start time in seconds since the epoch (midnight 1/1/1970) 4598 4599 # this is being used to seperate one number from a list. 4600 # what it is actually doing is sorting lists from numeric arrays. 4601 if type(times) is list or type(times) is ArrayType: 4602 number_of_times = len(times) 4603 times = ensure_numeric(times) 4604 if number_of_times == 0: 4605 starttime = 0 4606 else: 4607 starttime = times[0] 4608 times = times - starttime #Store relative times 4609 else: 4610 number_of_times = 0 4611 starttime = times 4612 #print "times",times 4613 #print "starttime", starttime 4614 outfile.starttime = starttime 4615 # dimension definitions 4616 outfile.createDimension('number_of_volumes', number_of_volumes) 4617 outfile.createDimension('number_of_vertices', 3) 4618 outfile.createDimension('numbers_in_range', 2) 4619 4620 if smoothing is True: 4621 outfile.createDimension('number_of_points', number_of_points) 4622 4623 # FIXME(Ole): This will cause sww files for paralle domains to 4624 # have ghost nodes stored (but not used by triangles). 4625 # To clean this up, we have to change get_vertex_values and 4626 # friends in quantity.py (but I can't be bothered right now) 4627 else: 4628 outfile.createDimension('number_of_points', 3*number_of_volumes) 4629 outfile.createDimension('number_of_timesteps', number_of_times) 4630 4631 # variable definitions 4632 outfile.createVariable('x', precision, ('number_of_points',)) 4633 outfile.createVariable('y', precision, ('number_of_points',)) 4634 outfile.createVariable('elevation', precision, ('number_of_points',)) 4635 q = 'elevation' 4636 outfile.createVariable(q+Write_sww.RANGE, precision, 4637 ('numbers_in_range',)) 4638 # The values are initally filled with large (10e+36) numbers. 4639 # I'm relying on this feature. Maybe I shouldn't? 4640 outfile.variables[q+Write_sww.RANGE][1] = \ 4641 -1*outfile.variables[q+Write_sww.RANGE][1] 4642 4643 #FIXME: Backwards compatibility 4644 outfile.createVariable('z', precision, ('number_of_points',)) 4645 ################################# 4646 4647 outfile.createVariable('volumes', Int, ('number_of_volumes', 4648 'number_of_vertices')) 4649 4650 outfile.createVariable('time', precision, 4651 ('number_of_timesteps',)) 4652 4653 for q in Write_sww.sww_quantities: 4654 outfile.createVariable(q, precision, 4655 ('number_of_timesteps', 4656 'number_of_points')) 4657 outfile.createVariable(q+Write_sww.RANGE, precision, 4658 ('numbers_in_range',)) 4659 # Initialising the max value to a very small number. 4660 # It assumes that netcdf initialises it to a very large number 4661 outfile.variables[q+Write_sww.RANGE][1] = \ 4662 -outfile.variables[q+Write_sww.RANGE][1] 4663 outfile.variables['time'][:] = times #Store time relative 4664 4665 if verbose: 4666 print '------------------------------------------------' 4667 print 'Statistics:' 4668 print ' t in [%f, %f], len(t) == %d'\ 4669 %(min(times.flat), max(times.flat), len(times.flat)) 4670 4671 4672 def triangulation(self,outfile, points_utm, volumes, 4673 elevation, zone=None, new_origin=None, 4674 points_georeference=None, verbose=False): 4675 """ 4676 4677 new_origin - qa georeference that the points can be set to. (Maybe 4678 do this before calling this function.) 4679 4680 points_utm - currently a list or array of the points in UTM. 4681 points_georeference - the georeference of the points_utm 4682 4683 How about passing new_origin and current_origin. 4684 If you get both, do a convertion from the old to the new. 4685 4686 If you only get new_origin, the points are absolute, convert to relative 4687 4688 if you only get the current_origin the points are relative, store 4689 as relative. 4690 4691 if you get no georefs create a new georef based on the minimums of 4692 points_utm. (Another option would be to default to absolute) 4693 4694 Yes, and this is done in another part of the code. 4695 Probably geospatial. 4696 4697 If you don't supply either geo_refs, then supply a zone. If not 4698 the default zone will be used. 4699 4700 4701 precon 4702 4703 header has been called. 4704 """ 4705 4706 number_of_points = len(points_utm) 4707 volumes = array(volumes) 4708 4709 # given the two geo_refs and the points, do the stuff 4710 # described in the method header 4711 # if this is needed else where, pull out as a function 4712 points_georeference = ensure_geo_reference(points_georeference) 4713 new_origin = ensure_geo_reference(new_origin) 4714 if new_origin is None and points_georeference is not None: 4715 points = points_utm 4716 geo_ref = points_georeference 4717 else: 4718 if new_origin is None: 4719 new_origin = Geo_reference(zone,min(points_utm[:,0]), 4720 min(points_utm[:,1])) 4721 points = new_origin.change_points_geo_ref(points_utm, 4722 points_georeference) 4723 geo_ref = new_origin 4724 4725 # At this stage I need a georef and points 4726 # the points are relative to the georef 4727 geo_ref.write_NetCDF(outfile) 4728 4729 # This will put the geo ref in the middle 4730 #geo_ref = Geo_reference(refzone,(max(x)+min(x))/2.0,(max(x)+min(y))/2.) 4731 x = points[:,0] 4732 y = points[:,1] 4733 z = outfile.variables['z'][:] 4734 4735 if verbose: 4736 print '------------------------------------------------' 4737 print 'More Statistics:' 4738 print ' Extent (/lon):' 4739 print ' x in [%f, %f], len(lat) == %d'\ 4740 %(min(x), max(x), 4741 len(x)) 4742 print ' y in [%f, %f], len(lon) == %d'\ 4743 %(min(y), max(y), 4744 len(y)) 4745 print ' z in [%f, %f], len(z) == %d'\ 4746 %(min(elevation), max(elevation), 4747 len(elevation)) 4748 print 'geo_ref: ',geo_ref 4749 print '------------------------------------------------' 4750 4751 #z = resize(bath_grid,outfile.variables['z'][:].shape) 4752 outfile.variables['x'][:] = points[:,0] #- geo_ref.get_xllcorner() 4753 outfile.variables['y'][:] = points[:,1] #- geo_ref.get_yllcorner() 4754 outfile.variables['z'][:] = elevation 4755 outfile.variables['elevation'][:] = elevation #FIXME HACK 4756 outfile.variables['volumes'][:] = volumes.astype(Int32) #On Opteron 64 4757 4758 q = 'elevation' 4759 # This updates the _range values 4760 outfile.variables[q+Write_sww.RANGE][0] = min(elevation) 4761 outfile.variables[q+Write_sww.RANGE][1] = max(elevation) 4762 4763 def quantities(self, outfile, precision=Float, 4764 slice_index=None, time=None, 4765 verbose=False, **quant): 4766 """ 4767 Write the quantity info. 4768 4769 **quant is extra keyword arguments passed in. These must be 4770 the sww quantities, currently; stage, xmomentum, ymomentum. 4771 4772 if the time array is already been built, use the slice_index 4773 to specify the index. 4774 4775 Otherwise, use time to increase the time dimension 4776 4777 Maybe make this general, but the viewer assumes these quantities, 4778 so maybe we don't want it general - unless the viewer is general 4779 4780 precon 4781 triangulation and 4782 header have been called. 4783 """ 4784 4785 if time is not None: 4786 file_time = outfile.variables['time'] 4787 slice_index = len(file_time) 4788 file_time[slice_index] = time 4789 4790 # write the conserved quantities from Doamin. 4791 # Typically stage, xmomentum, ymomentum 4792 # other quantities will be ignored, silently. 4793 for q in Write_sww.sww_quantities: 4794 if not quant.has_key(q): 4795 msg = 'SWW file can not write quantity %s' %q 4796 raise NewQuantity, msg 4797 else: 4798 q_values = quant[q] 4799 outfile.variables[q][slice_index] = q_values.astype(precision) 4800 4801 # This updates the _range values 4802 q_range = outfile.variables[q+Write_sww.RANGE][:] 4803 q_values_min = min(q_values) 4804 if q_values_min < q_range[0]: 4805 outfile.variables[q+Write_sww.RANGE][0] = q_values_min 4806 q_values_max = max(q_values) 4807 if q_values_max > q_range[1]: 4808 outfile.variables[q+Write_sww.RANGE][1] = q_values_max 4809 4810 def verbose_quantities(self, outfile): 4612 4811 print '------------------------------------------------' 4613 4812 print 'More Statistics:' 4614 print ' Extent (/lon):' 4615 print ' x in [%f, %f], len(lat) == %d'\ 4616 %(min(x), max(x), 4617 len(x)) 4618 print ' y in [%f, %f], len(lon) == %d'\ 4619 %(min(y), max(y), 4620 len(y)) 4621 print ' z in [%f, %f], len(z) == %d'\ 4622 %(min(elevation), max(elevation), 4623 len(elevation)) 4624 print 'geo_ref: ',geo_ref 4813 for q in Write_sww.sww_quantities: 4814 print ' %s in [%f, %f]' %(q 4815 ,outfile.variables[q+Write_sww.RANGE][0] 4816 ,outfile.variables[q+Write_sww.RANGE][1] 4817 ) 4625 4818 print '------------------------------------------------' 4626 4627 #z = resize(bath_grid,outfile.variables['z'][:].shape) 4628 outfile.variables['x'][:] = points_utm[:,0] - geo_ref.get_xllcorner() 4629 outfile.variables['y'][:] = points_utm[:,1] - geo_ref.get_yllcorner() 4630 outfile.variables['z'][:] = elevation 4631 outfile.variables['elevation'][:] = elevation #FIXME HACK 4632 outfile.variables['volumes'][:] = volumes.astype(Int32) #On Opteron 64 4633 4634 4635 def write_sww_time_slices(outfile, has, uas, vas, elevation, 4819 4820 def obsolete_write_sww_time_slices(outfile, has, uas, vas, elevation, 4636 4821 mean_stage=0, zscale=1, 4637 4822 verbose=False): … … 4651 4836 ymomentum[j] = -1*va*h # -1 since in mux files south is positive. 4652 4837 j += 1 4653 4654 def write_sww_time_slice(outfile, ha, ua, va, elevation, slice_index, 4655 mean_stage=0, zscale=1, 4656 verbose=False): 4657 #Time stepping 4658 stage = outfile.variables['stage'] 4659 xmomentum = outfile.variables['xmomentum'] 4660 ymomentum = outfile.variables['ymomentum'] 4661 4662 w = zscale*ha + mean_stage 4663 stage[slice_index] = w 4664 h = w - elevation 4665 xmomentum[slice_index] = ua*h 4666 ymomentum[slice_index] = -1*va*h # -1 since in mux files south is positive. 4667 4838 4668 4839 def urs2txt(basename_in, location_index=None): 4669 4840 """
Note: See TracChangeset
for help on using the changeset viewer.
