"""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 #Not yet implemented .nc: Native ferret NetCDF format FIXME: What else """ 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) 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 #Start time in seconds since the epoch (midnight 1/1/1970) fid.starttime = domain.starttime # 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 more general and rethink naming 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(filename, 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, ext = os.path.splitext(filename) #Get NetCDF infile = NetCDFFile(filename, 'r') #Open existing netcdf file for read if verbose: print 'Reading DEM from %s' %filename 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'] #Get output file xyaname = root + '.pts' if verbose: print 'Store to NetCDF file %s' %xyaname # NetCDF file definition outfile = NetCDFFile(xyaname, '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 = 56 #FIXME: Must be read from somewhere. Talk to Don. outfile.xllcorner = xllcorner #Easting of lower left corner outfile.yllcorner = yllcorner #Northing of lower left corner #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(filename, 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 to NetCDF format (.dem) mimcing the ASCII format closely. """ import os from Scientific.IO.NetCDF import NetCDFFile from Numeric import Float, array root, ext = os.path.splitext(filename) fid = open(filename) if verbose: print 'Reading DEM from %s' %filename lines = fid.readlines() fid.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 netcdfname = root + '.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 # 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(basefilename, verbose=False, minlat = None, maxlat =None, minlon = None, maxlon =None, mint = None, maxt = None, mean_stage = 0): """Convert 'Ferret' NetCDF format for wave propagation to sww format native to pyvolution. Specify only basefilename 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. """ import os from Scientific.IO.NetCDF import NetCDFFile from Numeric import Float, Int, searchsorted, zeros, array precision = Float #Get NetCDF data file_h = NetCDFFile(basefilename + '_ha.nc', 'r') #Wave amplitude (cm) file_u = NetCDFFile(basefilename + '_ua.nc', 'r') #Velocity (x) (cm/s) file_v = NetCDFFile(basefilename + '_va.nc', 'r') #Velocity (y) (cm/s) swwname = basefilename + '.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) latitudes = latitudes[kmin:kmax] longitudes = longitudes[lmin:lmax] times = times[jmin:jmax] 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_latitudes = latitudes.shape[0] number_of_longitudes = longitudes.shape[0] #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() 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'\ %(basefilename + '_ha.nc', basefilename + '_ua.nc', basefilename + '_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 l, lon in enumerate(longitudes): for k, lat in enumerate(latitudes): 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) origin = (min(x), min(y)) outfile.minx = origin[0] outfile.miny = origin[1] #print x - origin[0] #print y - origin[1] #print len(x) #print number_of_points outfile.variables['x'][:] = x - origin[0] outfile.variables['y'][:] = y - origin[1] outfile.variables['z'][:] = 0.0 outfile.variables['elevation'][:] = 0.0 outfile.variables['volumes'][:] = volumes outfile.variables['time'][:] = times #Time stepping stage = outfile.variables['stage'] xmomentum = outfile.variables['xmomentum'] ymomentum = outfile.variables['ymomentum'] for j in range(len(times)): i = 0 for l in range(number_of_longitudes): for k in range(number_of_latitudes): h = 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() #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()