"""Functions to store and retrieve data for the Shallow Water Wave equation. There are two kinds of data 1: Constant data: Vertex coordinates and field values. Stored once 2: Variable data: Conserved quantities. Stored once per timestep. All data is assumed to reside at vertex locations. Formats used within ANUGA: .sww: Netcdf format for storing model output .xya: ASCII format for storing arbitrary points and associated attributes .pts: NetCDF format for storing arbitrary points and associated attributes .asc: ASCII foramt of regular DEMs as output from ArcView .dem: NetCDF representation of regular DEM data .tsh: ASCII format for storing meshes and associated boundary and region info .nc: Native ferret NetCDF format A typical dataflow can be described as follows Manually created files: ASC, PRJ: Digital elevation models (gridded) TSH: Triangular meshes (e.g. dreated from pmesh) NC Model outputs for use as boundary conditions (e.g from MOST) AUTOMATICALLY CREATED FILES: ASC, PRJ -> DEM -> PTS: Conversion of DEM's to native pts file NC -> SWW: Conversion of MOST bundary files to boundary sww PTS + TSH -> TSH with elevation: Least squares fit TSH -> SWW: Conversion of TSH to sww viewable using Swollen TSH + Boundary SWW -> SWW: SImluation using pyvolution """ from Numeric import concatenate def make_filename(s): """Transform argument string into a suitable filename """ s = s.strip() s = s.replace(' ', '_') s = s.replace('(', '') s = s.replace(')', '') s = s.replace('__', '_') return s def check_dir(path, verbose=None): """Check that specified path exists. If path does not exist it will be created if possible USAGE: checkdir(path, verbose): ARGUMENTS: path -- Directory verbose -- Flag verbose output (default: None) RETURN VALUE: Verified path including trailing separator """ import os, sys import os.path if sys.platform in ['nt', 'dos', 'win32', 'what else?']: unix = 0 else: unix = 1 if path[-1] != os.sep: path = path + os.sep # Add separator for directories path = os.path.expanduser(path) # Expand ~ or ~user in pathname if not (os.access(path,os.R_OK and os.W_OK) or path == ''): try: exitcode=os.mkdir(path) # Change access rights if possible # if unix: exitcode=os.system('chmod 775 '+path) else: pass # FIXME: What about acces rights under Windows? if verbose: print 'MESSAGE: Directory', path, 'created.' except: print 'WARNING: Directory', path, 'could not be created.' if unix: path = '/tmp/' else: path = 'C:' print 'Using directory %s instead' %path return(path) def del_dir(path): """Recursively delete directory path and all its contents """ import os if os.path.isdir(path): for file in os.listdir(path): X = os.path.join(path, file) if os.path.isdir(X) and not os.path.islink(X): del_dir(X) else: try: os.remove(X) except: print "Could not remove file %s" %X os.rmdir(path) def create_filename(datadir, filename, format, size=None, time=None): import os #from config import data_dir FN = check_dir(datadir) + filename if size is not None: FN += '_size%d' %size if time is not None: FN += '_time%.2f' %time FN += '.' + format return FN def get_files(datadir, filename, format, size): """Get all file (names) with gven name, size and format """ import glob import os #from config import data_dir dir = check_dir(datadir) pattern = dir + os.sep + filename + '_size=%d*.%s' %(size, format) return glob.glob(pattern) #Generic class for storing output to e.g. visualisation or checkpointing class Data_format: """Generic interface to data formats """ def __init__(self, domain, extension, mode = 'w'): assert mode in ['r', 'w', 'a'], '''Mode %s must be either:''' %mode +\ ''' 'w' (write)'''+\ ''' 'r' (read)''' +\ ''' 'a' (append)''' #Create filename #self.filename = create_filename(domain.get_datadir(), #domain.get_name(), extension, len(domain)) self.filename = create_filename(domain.get_datadir(), domain.get_name(), extension) #print 'F', self.filename self.timestep = 0 self.number_of_volumes = len(domain) self.domain = domain #FIXME: Should we have a general set_precision function? #Class for storing output to e.g. visualisation class Data_format_sww(Data_format): """Interface to native NetCDF format (.sww) """ def __init__(self, domain, mode = 'w'): from Scientific.IO.NetCDF import NetCDFFile from Numeric import Int, Float, Float32 self.precision = Float32 #Use single precision Data_format.__init__(self, domain, 'sww', mode) # NetCDF file definition fid = NetCDFFile(self.filename, mode) if mode == 'w': #Create new file fid.institution = 'Geoscience Australia' fid.description = 'Output from pyvolution suitable for plotting' if domain.smooth: fid.smoothing = 'Yes' else: fid.smoothing = 'No' fid.order = domain.default_order #Reference point #Start time in seconds since the epoch (midnight 1/1/1970) fid.starttime = domain.starttime fid.xllcorner = domain.xllcorner fid.yllcorner = domain.yllcorner fid.zone = str(domain.zone) #FIXME: ? # dimension definitions fid.createDimension('number_of_volumes', self.number_of_volumes) fid.createDimension('number_of_vertices', 3) if domain.smooth is True: fid.createDimension('number_of_points', len(domain.vertexlist)) else: fid.createDimension('number_of_points', 3*self.number_of_volumes) fid.createDimension('number_of_timesteps', None) #extensible # variable definitions fid.createVariable('x', self.precision, ('number_of_points',)) fid.createVariable('y', self.precision, ('number_of_points',)) fid.createVariable('elevation', self.precision, ('number_of_points',)) #FIXME: Backwards compatibility fid.createVariable('z', self.precision, ('number_of_points',)) ################################# fid.createVariable('volumes', Int, ('number_of_volumes', 'number_of_vertices')) fid.createVariable('time', self.precision, ('number_of_timesteps',)) fid.createVariable('stage', self.precision, ('number_of_timesteps', 'number_of_points')) fid.createVariable('xmomentum', self.precision, ('number_of_timesteps', 'number_of_points')) fid.createVariable('ymomentum', self.precision, ('number_of_timesteps', 'number_of_points')) #Close fid.close() def store_connectivity(self): """Specialisation of store_connectivity for net CDF format Writes x,y,z coordinates of triangles constituting the bed elevation. """ from Scientific.IO.NetCDF import NetCDFFile from Numeric import concatenate, Int domain = self.domain #Get NetCDF fid = NetCDFFile(self.filename, 'a') #Open existing file for append # Get the variables x = fid.variables['x'] y = fid.variables['y'] z = fid.variables['elevation'] volumes = fid.variables['volumes'] # Get X, Y and bed elevation Z Q = domain.quantities['elevation'] X,Y,Z,V = Q.get_vertex_values(xy=True, precision = self.precision) x[:] = X.astype(self.precision) y[:] = Y.astype(self.precision) z[:] = Z.astype(self.precision) #FIXME: Backwards compatibility z = fid.variables['z'] z[:] = Z.astype(self.precision) ################################ volumes[:] = V.astype(volumes.typecode()) #Close fid.close() def store_timestep(self, names): """Store time and named quantities to file """ from Scientific.IO.NetCDF import NetCDFFile import types from time import sleep #Get NetCDF retries = 0 file_open = False while not file_open and retries < 10: try: fid = NetCDFFile(self.filename, 'a') #Open existing file except IOError: #This could happen if someone was reading the file. #In that case, wait a while and try again msg = 'Warning (store_timestep): File %s could not be opened'\ %self.filename msg += ' - trying again' print msg retries += 1 sleep(1) else: file_open = True if not file_open: msg = 'File %s could not be opened for append' %self.filename raise msg domain = self.domain # Get the variables time = fid.variables['time'] stage = fid.variables['stage'] xmomentum = fid.variables['xmomentum'] ymomentum = fid.variables['ymomentum'] i = len(time) #Store time time[i] = self.domain.time if type(names) not in [types.ListType, types.TupleType]: names = [names] for name in names: # Get quantity Q = domain.quantities[name] A,V = Q.get_vertex_values(xy=False, precision = self.precision) #FIXME: Make this general (see below) if name == 'stage': stage[i,:] = A.astype(self.precision) elif name == 'xmomentum': xmomentum[i,:] = A.astype(self.precision) elif name == 'ymomentum': ymomentum[i,:] = A.astype(self.precision) #As in.... #eval( name + '[i,:] = A.astype(self.precision)' ) #FIXME: But we need a UNIT test for that before refactoring #Flush and close fid.sync() fid.close() #Class for handling checkpoints data class Data_format_cpt(Data_format): """Interface to native NetCDF format (.cpt) """ def __init__(self, domain, mode = 'w'): from Scientific.IO.NetCDF import NetCDFFile from Numeric import Int, Float, Float self.precision = Float #Use full precision Data_format.__init__(self, domain, 'sww', mode) # NetCDF file definition fid = NetCDFFile(self.filename, mode) if mode == 'w': #Create new file fid.institution = 'Geoscience Australia' fid.description = 'Checkpoint data' #fid.smooth = domain.smooth fid.order = domain.default_order # dimension definitions fid.createDimension('number_of_volumes', self.number_of_volumes) fid.createDimension('number_of_vertices', 3) #Store info at all vertices (no smoothing) fid.createDimension('number_of_points', 3*self.number_of_volumes) fid.createDimension('number_of_timesteps', None) #extensible # variable definitions #Mesh fid.createVariable('x', self.precision, ('number_of_points',)) fid.createVariable('y', self.precision, ('number_of_points',)) fid.createVariable('volumes', Int, ('number_of_volumes', 'number_of_vertices')) fid.createVariable('time', self.precision, ('number_of_timesteps',)) #Allocate space for all quantities for name in domain.quantities.keys(): fid.createVariable(name, self.precision, ('number_of_timesteps', 'number_of_points')) #Close fid.close() def store_checkpoint(self): """ Write x,y coordinates of triangles. Write connectivity ( constituting the bed elevation. """ from Scientific.IO.NetCDF import NetCDFFile from Numeric import concatenate domain = self.domain #Get NetCDF fid = NetCDFFile(self.filename, 'a') #Open existing file for append # Get the variables x = fid.variables['x'] y = fid.variables['y'] volumes = fid.variables['volumes'] # Get X, Y and bed elevation Z Q = domain.quantities['elevation'] X,Y,Z,V = Q.get_vertex_values(xy=True, precision = self.precision) x[:] = X.astype(self.precision) y[:] = Y.astype(self.precision) z[:] = Z.astype(self.precision) volumes[:] = V #Close fid.close() def store_timestep(self, name): """Store time and named quantity to file """ from Scientific.IO.NetCDF import NetCDFFile from time import sleep #Get NetCDF retries = 0 file_open = False while not file_open and retries < 10: try: fid = NetCDFFile(self.filename, 'a') #Open existing file except IOError: #This could happen if someone was reading the file. #In that case, wait a while and try again msg = 'Warning (store_timestep): File %s could not be opened'\ %self.filename msg += ' - trying again' print msg retries += 1 sleep(1) else: file_open = True if not file_open: msg = 'File %s could not be opened for append' %self.filename raise msg domain = self.domain # Get the variables time = fid.variables['time'] stage = fid.variables['stage'] i = len(time) #Store stage time[i] = self.domain.time # Get quantity Q = domain.quantities[name] A,V = Q.get_vertex_values(xy=False, precision = self.precision) stage[i,:] = A.astype(self.precision) #Flush and close fid.sync() fid.close() #Function for storing xya output #FIXME Not done yet for this version class Data_format_xya(Data_format): """Generic interface to data formats """ def __init__(self, domain, mode = 'w'): from Scientific.IO.NetCDF import NetCDFFile from Numeric import Int, Float, Float32 self.precision = Float32 #Use single precision Data_format.__init__(self, domain, 'xya', mode) #FIXME -This is the old xya format def store_all(self): """Specialisation of store all for xya format Writes x,y,z coordinates of triangles constituting the bed elevation. """ from Numeric import concatenate domain = self.domain fd = open(self.filename, 'w') if domain.smooth is True: number_of_points = len(domain.vertexlist) else: number_of_points = 3*self.number_of_volumes numVertAttrib = 3 #Three attributes is what is assumed by the xya format fd.write(str(number_of_points) + " " + str(numVertAttrib) +\ " # [attributes]" + "\n") # Get X, Y, bed elevation and friction (index=0,1) X,Y,A,V = domain.get_vertex_values(xy=True, value_array='field_values', indices = (0,1), precision = self.precision) bed_eles = A[:,0] fricts = A[:,1] # Get stage (index=0) B,V = domain.get_vertex_values(xy=False, value_array='conserved_quantities', indices = (0,), precision = self.precision) stages = B[:,0] # [attributes] for x, y, bed_ele, stage, frict in map(None, X, Y, bed_eles, stages, fricts): s = '%.6f %.6f %.6f %.6f %.6f\n' %(x, y, bed_ele, stage, frict) fd.write(s) #close fd.close() def store_timestep(self, t, V0, V1, V2): """Store time, water heights (and momentums) to file """ pass #Auxiliary def write_obj(filename,x,y,z): """Store x,y,z vectors into filename (obj format) Vectors are assumed to have dimension (M,3) where M corresponds to the number elements. triangles are assumed to be disconnected The three numbers in each vector correspond to three vertices, e.g. the x coordinate of vertex 1 of element i is in x[i,1] """ #print 'Writing obj to %s' % filename import os.path root, ext = os.path.splitext(filename) if ext == '.obj': FN = filename else: FN = filename + '.obj' outfile = open(FN, 'wb') outfile.write("# Triangulation as an obj file\n") M, N = x.shape assert N==3 #Assuming three vertices per element for i in range(M): for j in range(N): outfile.write("v %f %f %f\n" % (x[i,j],y[i,j],z[i,j])) for i in range(M): base = i*N outfile.write("f %d %d %d\n" % (base+1,base+2,base+3)) outfile.close() #Conversion routines def sww2obj(basefilename, size): """Convert netcdf based data output to obj """ from Scientific.IO.NetCDF import NetCDFFile from Numeric import Float, zeros #Get NetCDF FN = create_filename('.', basefilename, 'sww', size) print 'Reading from ', FN fid = NetCDFFile(FN, 'r') #Open existing file for read # Get the variables x = fid.variables['x'] y = fid.variables['y'] z = fid.variables['elevation'] time = fid.variables['time'] stage = fid.variables['stage'] M = size #Number of lines xx = zeros((M,3), Float) yy = zeros((M,3), Float) zz = zeros((M,3), Float) for i in range(M): for j in range(3): xx[i,j] = x[i+j*M] yy[i,j] = y[i+j*M] zz[i,j] = z[i+j*M] #Write obj for bathymetry FN = create_filename('.', basefilename, 'obj', size) write_obj(FN,xx,yy,zz) #Now read all the data with variable information, combine with #x,y info and store as obj for k in range(len(time)): t = time[k] print 'Processing timestep %f' %t for i in range(M): for j in range(3): zz[i,j] = stage[k,i+j*M] #Write obj for variable data #FN = create_filename(basefilename, 'obj', size, time=t) FN = create_filename('.', basefilename[:5], 'obj', size, time=t) write_obj(FN,xx,yy,zz) def dat2obj(basefilename): """Convert line based data output to obj FIXME: Obsolete? """ import glob, os from config import data_dir #Get bathymetry and x,y's lines = open(data_dir+os.sep+basefilename+'_geometry.dat', 'r').readlines() from Numeric import zeros, Float M = len(lines) #Number of lines x = zeros((M,3), Float) y = zeros((M,3), Float) z = zeros((M,3), Float) ##i = 0 for i, line in enumerate(lines): tokens = line.split() values = map(float,tokens) for j in range(3): x[i,j] = values[j*3] y[i,j] = values[j*3+1] z[i,j] = values[j*3+2] ##i += 1 #Write obj for bathymetry write_obj(data_dir+os.sep+basefilename+'_geometry',x,y,z) #Now read all the data files with variable information, combine with #x,y info #and store as obj files = glob.glob(data_dir+os.sep+basefilename+'*.dat') for filename in files: print 'Processing %s' % filename lines = open(data_dir+os.sep+filename,'r').readlines() assert len(lines) == M root, ext = os.path.splitext(filename) #Get time from filename i0 = filename.find('_time=') if i0 == -1: #Skip bathymetry file continue i0 += 6 #Position where time starts i1 = filename.find('.dat') if i1 > i0: t = float(filename[i0:i1]) else: raise 'Hmmmm' ##i = 0 for i, line in enumerate(lines): tokens = line.split() values = map(float,tokens) for j in range(3): z[i,j] = values[j] ##i += 1 #Write obj for variable data write_obj(data_dir+os.sep+basefilename+'_time=%.4f' %t,x,y,z) def filter_netcdf(filename1, filename2, first=0, last=None, step = 1): """Read netcdf filename1, pick timesteps first:step:last and save to nettcdf file filename2 """ from Scientific.IO.NetCDF import NetCDFFile #Get NetCDF infile = NetCDFFile(filename1, 'r') #Open existing file for read outfile = NetCDFFile(filename2, 'w') #Open new file #Copy dimensions for d in infile.dimensions: outfile.createDimension(d, infile.dimensions[d]) for name in infile.variables: var = infile.variables[name] outfile.createVariable(name, var.typecode(), var.dimensions) #Copy the static variables for name in infile.variables: if name == 'time' or name == 'stage': pass else: #Copy outfile.variables[name][:] = infile.variables[name][:] #Copy selected timesteps time = infile.variables['time'] stage = infile.variables['stage'] newtime = outfile.variables['time'] newstage = outfile.variables['stage'] if last is None: last = len(time) selection = range(first, last, step) for i, j in enumerate(selection): print 'Copying timestep %d of %d (%f)' %(j, last-first, time[j]) newtime[i] = time[j] newstage[i,:] = stage[j,:] #Close infile.close() outfile.close() #Get data objects def get_dataobject(domain, mode='w'): """Return instance of class of given format using filename """ cls = eval('Data_format_%s' %domain.format) return cls(domain, mode) def dem2pts(basename_in, basename_out=None, verbose=False): """Read Digitial Elevation model from the following NetCDF format (.dem) Example: ncols 3121 nrows 1800 xllcorner 722000 yllcorner 5893000 cellsize 25 NODATA_value -9999 138.3698 137.4194 136.5062 135.5558 .......... Convert to NetCDF pts format which is points: (Nx2) Float array elevation: N Float array """ import os from Scientific.IO.NetCDF import NetCDFFile from Numeric import Float, arrayrange, concatenate root = basename_in #Get NetCDF infile = NetCDFFile(root + '.dem', 'r') #Open existing netcdf file for read if verbose: print 'Reading DEM from %s' %(root + '.dem') ncols = infile.ncols[0] nrows = infile.nrows[0] xllcorner = infile.xllcorner[0] #Easting of lower left corner yllcorner = infile.yllcorner[0] #Northing of lower left corner cellsize = infile.cellsize[0] NODATA_value = infile.NODATA_value[0] dem_elevation = infile.variables['elevation'] zone = infile.zone[0] false_easting = infile.false_easting[0] false_northing = infile.false_northing[0] #Text strings projection = infile.projection datum = infile.datum units = infile.units #Get output file if basename_out == None: ptsname = root + '.pts' else: ptsname = basename_out + '.pts' if verbose: print 'Store to NetCDF file %s' %ptsname # NetCDF file definition outfile = NetCDFFile(ptsname, 'w') #Create new file outfile.institution = 'Geoscience Australia' outfile.description = 'NetCDF pts format for compact and portable storage ' +\ 'of spatial point data' #Georeferencing outfile.zone = zone outfile.xllcorner = xllcorner #Easting of lower left corner outfile.yllcorner = yllcorner #Northing of lower left corner outfile.false_easting = false_easting outfile.false_northing =false_northing outfile.projection = projection outfile.datum = datum outfile.units = units #Grid info (FIXME: probably not going to be used, but heck) outfile.ncols = ncols outfile.nrows = nrows # dimension definitions outfile.createDimension('number_of_points', nrows*ncols) outfile.createDimension('number_of_dimensions', 2) #This is 2d data # variable definitions outfile.createVariable('points', Float, ('number_of_points', 'number_of_dimensions')) outfile.createVariable('elevation', Float, ('number_of_points',)) # Get handles to the variables points = outfile.variables['points'] elevation = outfile.variables['elevation'] #Store data #FIXME: Could perhaps be faster using array operations for i in range(nrows): if verbose: print 'Processing row %d of %d' %(i, nrows) y = (nrows-i)*cellsize for j in range(ncols): index = i*ncols + j x = j*cellsize points[index, :] = [x,y] elevation[index] = dem_elevation[i, j] infile.close() outfile.close() def convert_dem_from_ascii2netcdf(basename_in, basename_out = None, verbose=False): """Read Digitial Elevation model from the following ASCII format (.asc) Example: ncols 3121 nrows 1800 xllcorner 722000 yllcorner 5893000 cellsize 25 NODATA_value -9999 138.3698 137.4194 136.5062 135.5558 .......... Convert basename_in + '.asc' to NetCDF format (.dem) mimicking the ASCII format closely. An accompanying file with same basename_in but extension .prj must exist and is used to fix the UTM zone, datum, false northings and eastings. The prj format is assumed to be as Projection UTM Zone 56 Datum WGS84 Zunits NO Units METERS Spheroid WGS84 Xshift 0.0000000000 Yshift 10000000.0000000000 Parameters """ import os from Scientific.IO.NetCDF import NetCDFFile from Numeric import Float, array #root, ext = os.path.splitext(basename_in) root = basename_in ########################################### # Read Meta data if verbose: print 'Reading METADATA from %s' %root + '.prj' metadatafile = open(root + '.prj') metalines = metadatafile.readlines() metadatafile.close() L = metalines[0].strip().split() assert L[0].strip().lower() == 'projection' projection = L[1].strip() #TEXT L = metalines[1].strip().split() assert L[0].strip().lower() == 'zone' zone = int(L[1].strip()) L = metalines[2].strip().split() assert L[0].strip().lower() == 'datum' datum = L[1].strip() #TEXT L = metalines[3].strip().split() assert L[0].strip().lower() == 'zunits' #IGNORE zunits = L[1].strip() #TEXT L = metalines[4].strip().split() assert L[0].strip().lower() == 'units' units = L[1].strip() #TEXT L = metalines[5].strip().split() assert L[0].strip().lower() == 'spheroid' #IGNORE spheroid = L[1].strip() #TEXT L = metalines[6].strip().split() assert L[0].strip().lower() == 'xshift' false_easting = float(L[1].strip()) L = metalines[7].strip().split() assert L[0].strip().lower() == 'yshift' false_northing = float(L[1].strip()) #print false_easting, false_northing, zone, datum ########################################### #Read DEM data datafile = open(basename_in + '.asc') if verbose: print 'Reading DEM from %s' %(basename_in + '.asc') lines = datafile.readlines() datafile.close() if verbose: print 'Got', len(lines), ' lines' ncols = int(lines[0].split()[1].strip()) nrows = int(lines[1].split()[1].strip()) xllcorner = float(lines[2].split()[1].strip()) yllcorner = float(lines[3].split()[1].strip()) cellsize = float(lines[4].split()[1].strip()) NODATA_value = int(lines[5].split()[1].strip()) assert len(lines) == nrows + 6 ########################################## if basename_out == None: netcdfname = root + '.dem' else: netcdfname = basename_out + '.dem' if verbose: print 'Store to NetCDF file %s' %netcdfname # NetCDF file definition fid = NetCDFFile(netcdfname, 'w') #Create new file fid.institution = 'Geoscience Australia' fid.description = 'NetCDF DEM format for compact and portable storage ' +\ 'of spatial point data' fid.ncols = ncols fid.nrows = nrows fid.xllcorner = xllcorner fid.yllcorner = yllcorner fid.cellsize = cellsize fid.NODATA_value = NODATA_value fid.zone = zone fid.false_easting = false_easting fid.false_northing = false_northing fid.projection = projection fid.datum = datum fid.units = units # dimension definitions fid.createDimension('number_of_rows', nrows) fid.createDimension('number_of_columns', ncols) # variable definitions fid.createVariable('elevation', Float, ('number_of_rows', 'number_of_columns')) # Get handles to the variables elevation = fid.variables['elevation'] #Store data for i, line in enumerate(lines[6:]): fields = line.split() if verbose: print 'Processing row %d of %d' %(i, nrows) elevation[i, :] = array([float(x) for x in fields]) fid.close() def ferret2sww(basename_in, basename_out = None, verbose = False, minlat = None, maxlat = None, minlon = None, maxlon = None, mint = None, maxt = None, mean_stage = 0, origin = None, zscale = 1): """Convert 'Ferret' NetCDF format for wave propagation to sww format native to pyvolution. Specify only basename_in and read files of the form basefilename_ha.nc, basefilename_ua.nc, basefilename_va.nc containing relative height, x-velocity and y-velocity, respectively. Also convert latitude and longitude to UTM. All coordinates are assumed to be given in the GDA94 datum. min's and max's: If omitted - full extend is used. To include a value min may equal it, while max must exceed it. Lat and lon are assuemd to be in decimal degrees origin is a 3-tuple with geo referenced UTM coordinates (zone, easting, northing) nc format has values organised as HA[TIME, LATITUDE, LONGITUDE] which means that longitude is the fastest varying dimension (row major order, so to speak) ferret2sww uses grid points as vertices in a triangular grid counting vertices from lower left corner upwards, then right """ import os from Scientific.IO.NetCDF import NetCDFFile from Numeric import Float, Int, Int32, searchsorted, zeros, array precision = Float #Get NetCDF data if verbose: print 'Reading files %s_*.nc' %basename_in file_h = NetCDFFile(basename_in + '_ha.nc', 'r') #Wave amplitude (cm) file_u = NetCDFFile(basename_in + '_ua.nc', 'r') #Velocity (x) (cm/s) file_v = NetCDFFile(basename_in + '_va.nc', 'r') #Velocity (y) (cm/s) if basename_out is None: swwname = basename_in + '.sww' else: swwname = basename_out + '.sww' times = file_h.variables['TIME'] latitudes = file_h.variables['LAT'] longitudes = file_h.variables['LON'] if mint == None: jmin = 0 else: jmin = searchsorted(times, mint) if maxt == None: jmax=len(times) else: jmax = searchsorted(times, maxt) if minlat == None: kmin=0 else: kmin = searchsorted(latitudes, minlat) if maxlat == None: kmax = len(latitudes) else: kmax = searchsorted(latitudes, maxlat) if minlon == None: lmin=0 else: lmin = searchsorted(longitudes, minlon) if maxlon == None: lmax = len(longitudes) else: lmax = searchsorted(longitudes, maxlon) times = times[jmin:jmax] latitudes = latitudes[kmin:kmax] longitudes = longitudes[lmin:lmax] if verbose: print 'cropping' amplitudes = file_h.variables['HA'][jmin:jmax, kmin:kmax, lmin:lmax] xspeed = file_u.variables['UA'][jmin:jmax, kmin:kmax, lmin:lmax] yspeed = file_v.variables['VA'][jmin:jmax, kmin:kmax, lmin:lmax] number_of_times = times.shape[0] number_of_latitudes = latitudes.shape[0] number_of_longitudes = longitudes.shape[0] assert amplitudes.shape[0] == number_of_times assert amplitudes.shape[1] == number_of_latitudes assert amplitudes.shape[2] == number_of_longitudes #print times #print latitudes #print longitudes #print 'MIN', min(min(min(amplitudes))) #print 'MAX', max(max(max(amplitudes))) #print number_of_latitudes, number_of_longitudes number_of_points = number_of_latitudes*number_of_longitudes number_of_volumes = (number_of_latitudes-1)*(number_of_longitudes-1)*2 #print file_h.dimensions.keys() #print file_h.variables.keys() file_h.close() file_u.close() file_v.close() if verbose: print 'Store to SWW file %s' %swwname # NetCDF file definition outfile = NetCDFFile(swwname, 'w') #Create new file outfile.institution = 'Geoscience Australia' outfile.description = 'Converted from Ferret files: %s, %s, %s'\ %(basename_in + '_ha.nc', basename_in + '_ua.nc', basename_in + '_va.nc') #For sww compatibility outfile.smoothing = 'Yes' outfile.order = 1 #Start time in seconds since the epoch (midnight 1/1/1970) outfile.starttime = times[0] # dimension definitions outfile.createDimension('number_of_volumes', number_of_volumes) outfile.createDimension('number_of_vertices', 3) outfile.createDimension('number_of_points', number_of_points) #outfile.createDimension('number_of_timesteps', len(times)) outfile.createDimension('number_of_timesteps', len(times)) # variable definitions outfile.createVariable('x', precision, ('number_of_points',)) outfile.createVariable('y', precision, ('number_of_points',)) outfile.createVariable('elevation', precision, ('number_of_points',)) #FIXME: Backwards compatibility outfile.createVariable('z', precision, ('number_of_points',)) ################################# outfile.createVariable('volumes', Int, ('number_of_volumes', 'number_of_vertices')) outfile.createVariable('time', precision, ('number_of_timesteps',)) outfile.createVariable('stage', precision, ('number_of_timesteps', 'number_of_points')) outfile.createVariable('xmomentum', precision, ('number_of_timesteps', 'number_of_points')) outfile.createVariable('ymomentum', precision, ('number_of_timesteps', 'number_of_points')) #Store from coordinate_transforms.redfearn import redfearn x = zeros(number_of_points, Float) #Easting y = zeros(number_of_points, Float) #Northing #volumes = zeros(number_of_volumes, Int) i = 0 #Check zone boundaries refzone, _, _ = redfearn(latitudes[0],longitudes[0]) vertices = {} for k, lat in enumerate(latitudes): for l, lon in enumerate(longitudes): vertices[l,k] = i zone, easting, northing = redfearn(lat,lon) msg = 'Zone boundary crossed at longitude =', lon assert zone == refzone, msg #print '%7.2f %7.2f %8.2f %8.2f' %(lon, lat, easting, northing) x[i] = easting y[i] = northing i += 1 #Construct 2 triangles per 'rectangular' element volumes = [] for l in range(number_of_longitudes-1): for k in range(number_of_latitudes-1): v1 = vertices[l,k+1] v2 = vertices[l,k] v3 = vertices[l+1,k+1] v4 = vertices[l+1,k] volumes.append([v1,v2,v3]) #Upper element volumes.append([v4,v3,v2]) #Lower element volumes = array(volumes) if origin == None: zone = refzone xllcorner = min(x) yllcorner = min(y) else: zone = origin[0] xllcorner = origin[1] yllcorner = origin[2] outfile.xllcorner = xllcorner outfile.yllcorner = yllcorner outfile.zone = zone outfile.variables['x'][:] = x - xllcorner outfile.variables['y'][:] = y - yllcorner outfile.variables['z'][:] = 0.0 outfile.variables['elevation'][:] = 0.0 outfile.variables['time'][:] = times outfile.variables['volumes'][:] = volumes.astype(Int32) #On Opteron 64 #Time stepping stage = outfile.variables['stage'] xmomentum = outfile.variables['xmomentum'] ymomentum = outfile.variables['ymomentum'] for j in range(len(times)): i = 0 for k in range(number_of_latitudes): for l in range(number_of_longitudes): h = zscale*amplitudes[j,k,l]/100 + mean_stage stage[j,i] = h xmomentum[j,i] = xspeed[j,k,l]/100*h ymomentum[j,i] = yspeed[j,k,l]/100*h i += 1 outfile.close() def extent_sww(file_name): """ Read in an sww file. Input; file_name - the sww file Output; z - Vector of bed elevation volumes - Array. Each row has 3 values, representing the vertices that define the volume time - Vector of the times where there is stage information stage - array with respect to time and vertices (x,y) """ from Scientific.IO.NetCDF import NetCDFFile #Check contents #Get NetCDF fid = NetCDFFile(file_name, 'r') # Get the variables x = fid.variables['x'][:] y = fid.variables['y'][:] stage = fid.variables['stage'][:] #print "stage",stage #print "stage.shap",stage.shape #print "min(stage.flat), mpythonax(stage.flat)",min(stage.flat), max(stage.flat) #print "min(stage)",min(stage) fid.close() return [min(x),max(x),min(y),max(y),min(stage.flat),max(stage.flat)] #OBSOLETE STUFF #Native checkpoint format. #Information needed to recreate a state is preserved #FIXME: Rethink and maybe use netcdf format def cpt_variable_writer(filename, t, v0, v1, v2): """Store all conserved quantities to file """ M, N = v0.shape FN = create_filename(filename, 'cpt', M, t) #print 'Writing to %s' %FN fid = open(FN, 'w') for i in range(M): for j in range(N): fid.write('%.16e ' %v0[i,j]) for j in range(N): fid.write('%.16e ' %v1[i,j]) for j in range(N): fid.write('%.16e ' %v2[i,j]) fid.write('\n') fid.close() def cpt_variable_reader(filename, t, v0, v1, v2): """Store all conserved quantities to file """ M, N = v0.shape FN = create_filename(filename, 'cpt', M, t) #print 'Reading from %s' %FN fid = open(FN) for i in range(M): values = fid.readline().split() #Get one line for j in range(N): v0[i,j] = float(values[j]) v1[i,j] = float(values[3+j]) v2[i,j] = float(values[6+j]) fid.close() def cpt_constant_writer(filename, X0, X1, X2, v0, v1, v2): """Writes x,y,z,z,z coordinates of triangles constituting the bed elevation. Not in use pt """ M, N = v0.shape print X0 import sys; sys.exit() FN = create_filename(filename, 'cpt', M) print 'Writing to %s' %FN fid = open(FN, 'w') for i in range(M): for j in range(2): fid.write('%.16e ' %X0[i,j]) #x, y for j in range(N): fid.write('%.16e ' %v0[i,j]) #z,z,z, for j in range(2): fid.write('%.16e ' %X1[i,j]) #x, y for j in range(N): fid.write('%.16e ' %v1[i,j]) for j in range(2): fid.write('%.16e ' %X2[i,j]) #x, y for j in range(N): fid.write('%.16e ' %v2[i,j]) fid.write('\n') fid.close() #Function for storing out to e.g. visualisation #FIXME: Do we want this? #FIXME: Not done yet for this version def dat_constant_writer(filename, X0, X1, X2, v0, v1, v2): """Writes x,y,z coordinates of triangles constituting the bed elevation. """ M, N = v0.shape FN = create_filename(filename, 'dat', M) #print 'Writing to %s' %FN fid = open(FN, 'w') for i in range(M): for j in range(2): fid.write('%f ' %X0[i,j]) #x, y fid.write('%f ' %v0[i,0]) #z for j in range(2): fid.write('%f ' %X1[i,j]) #x, y fid.write('%f ' %v1[i,0]) #z for j in range(2): fid.write('%f ' %X2[i,j]) #x, y fid.write('%f ' %v2[i,0]) #z fid.write('\n') fid.close() def dat_variable_writer(filename, t, v0, v1, v2): """Store water height to file """ M, N = v0.shape FN = create_filename(filename, 'dat', M, t) #print 'Writing to %s' %FN fid = open(FN, 'w') for i in range(M): fid.write('%.4f ' %v0[i,0]) fid.write('%.4f ' %v1[i,0]) fid.write('%.4f ' %v2[i,0]) fid.write('\n') fid.close()