"""datamanager.py - input output for AnuGA This module takes care of reading and writing datafiles such as topograhies, model output, etc Formats used within AnuGA: .sww: Netcdf format for storing model output f(t,x,y) .tms: Netcdf format for storing time series f(t) .csv: ASCII format for storing arbitrary points and associated attributes .pts: NetCDF format for storing arbitrary points and associated attributes .asc: ASCII format of regular DEMs as output from ArcView .prj: Associated ArcView file giving more meta data for asc format .ers: ERMapper header format of regular DEMs for ArcView .dem: NetCDF representation of regular DEM data .tsh: ASCII format for storing meshes and associated boundary and region info .msh: NetCDF format for storing meshes and associated boundary and region info .nc: Native ferret NetCDF format .geo: Houdinis ascii geometry format (?) A typical dataflow can be described as follows Manually created files: ASC, PRJ: Digital elevation models (gridded) TSH: Triangular meshes (e.g. created from anuga.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 abstract_2d_finite_volumes """ import exceptions class TitleValueError(exceptions.Exception): pass class DataMissingValuesError(exceptions.Exception): pass class DataFileNotOpenError(exceptions.Exception): pass class DataTimeError(exceptions.Exception): pass class DataDomainError(exceptions.Exception): pass class NewQuantity(exceptions.Exception): pass import csv import os, sys import shutil from struct import unpack import array as p_array #import time, os from os import sep, path, remove, mkdir, access, F_OK, W_OK, getcwd from Numeric import concatenate, array, Float, Int, Int32, resize, sometrue, \ searchsorted, zeros, allclose, around, reshape, transpose, sort, \ NewAxis, ArrayType, compress, take, arange, argmax, alltrue,shape,Float import string from Scientific.IO.NetCDF import NetCDFFile #from shutil import copy from os.path import exists, basename, join from os import getcwd from anuga.coordinate_transforms.redfearn import redfearn, \ convert_from_latlon_to_utm from anuga.coordinate_transforms.geo_reference import Geo_reference, \ write_NetCDF_georeference, ensure_geo_reference from anuga.geospatial_data.geospatial_data import Geospatial_data,\ ensure_absolute from anuga.config import minimum_storable_height as default_minimum_storable_height from anuga.config import max_float from anuga.utilities.numerical_tools import ensure_numeric, mean from anuga.caching.caching import myhash from anuga.utilities.anuga_exceptions import ANUGAError from anuga.shallow_water import Domain from anuga.abstract_2d_finite_volumes.pmesh2domain import \ pmesh_to_domain_instance from anuga.abstract_2d_finite_volumes.util import get_revision_number, \ remove_lone_verts, sww2timeseries, get_centroid_values from anuga.load_mesh.loadASCII import export_mesh_file # formula mappings quantity_formula = {'momentum':'(xmomentum**2 + ymomentum**2)**0.5', 'depth':'stage-elevation', 'speed': \ '(xmomentum**2 + ymomentum**2)**0.5/(stage-elevation+1.e-6/(stage-elevation))'} def make_filename(s): """Transform argument string into a Sexsuitable 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.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) # ANOTHER OPTION, IF NEED IN THE FUTURE, Nick B 7/2007 def rmgeneric(path, __func__,verbose=False): ERROR_STR= """Error removing %(path)s, %(error)s """ try: __func__(path) if verbose: print 'Removed ', path except OSError, (errno, strerror): print ERROR_STR % {'path' : path, 'error': strerror } def removeall(path,verbose=False): if not os.path.isdir(path): return files=os.listdir(path) for x in files: fullpath=os.path.join(path, x) if os.path.isfile(fullpath): f=os.remove rmgeneric(fullpath, f) elif os.path.isdir(fullpath): removeall(fullpath) f=os.rmdir rmgeneric(fullpath, f,verbose) def create_filename(datadir, filename, format, size=None, time=None): import os #from anuga.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 given name, size and format """ import glob import os #from anuga.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) #print 'F', self.filename self.timestep = 0 self.domain = domain # Exclude ghosts in case this is a parallel domain self.number_of_nodes = domain.number_of_full_nodes self.number_of_volumes = domain.number_of_full_triangles #self.number_of_volumes = len(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) for storing model output 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. """ def __init__(self, domain, mode = 'w',\ max_size = 2000000000, recursion = False): from Scientific.IO.NetCDF import NetCDFFile from Numeric import Int, Float, Float32 self.precision = Float32 #Use single precision for quantities if hasattr(domain, 'max_size'): self.max_size = domain.max_size #file size max is 2Gig else: self.max_size = max_size self.recursion = recursion self.mode = mode Data_format.__init__(self, domain, 'sww', mode) if hasattr(domain, 'minimum_storable_height'): self.minimum_storable_height = domain.minimum_storable_height else: self.minimum_storable_height = default_minimum_storable_height # NetCDF file definition fid = NetCDFFile(self.filename, mode) if mode == 'w': description = 'Output from anuga.abstract_2d_finite_volumes suitable for plotting' self.writer = Write_sww() self.writer.store_header(fid, domain.starttime, self.number_of_volumes, self.domain.number_of_full_nodes, description=description, smoothing=domain.smooth, order=domain.default_order) # Extra optional information if hasattr(domain, 'texture'): fid.texture = domain.texture if domain.quantities_to_be_monitored is not None: fid.createDimension('singleton', 1) fid.createDimension('two', 2) poly = domain.monitor_polygon if poly is not None: N = len(poly) fid.createDimension('polygon_length', N) fid.createVariable('extrema.polygon', self.precision, ('polygon_length', 'two')) fid.variables['extrema.polygon'][:] = poly interval = domain.monitor_time_interval if interval is not None: fid.createVariable('extrema.time_interval', self.precision, ('two',)) fid.variables['extrema.time_interval'][:] = interval for q in domain.quantities_to_be_monitored: #print 'doing', q fid.createVariable(q+'.extrema', self.precision, ('numbers_in_range',)) fid.createVariable(q+'.min_location', self.precision, ('numbers_in_range',)) fid.createVariable(q+'.max_location', self.precision, ('numbers_in_range',)) fid.createVariable(q+'.min_time', self.precision, ('singleton',)) fid.createVariable(q+'.max_time', self.precision, ('singleton',)) 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) # points = concatenate( (X[:,NewAxis],Y[:,NewAxis]), axis=1 ) self.writer.store_triangulation(fid, points, V.astype(volumes.typecode()), Z, points_georeference= \ domain.geo_reference) # 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 from os import stat from Numeric import choose # 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 step %s again' %self.domain.time 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 DataFileNotOpenError, msg # Check to see if the file is already too big: time = fid.variables['time'] i = len(time)+1 file_size = stat(self.filename)[6] file_size_increase = file_size/i if file_size + file_size_increase > self.max_size*(2**self.recursion): # In order to get the file name and start time correct, # I change the domain.filename and domain.starttime. # This is the only way to do this without changing # other modules (I think). # Write a filename addon that won't break swollens reader # (10.sww is bad) filename_ext = '_time_%s'%self.domain.time filename_ext = filename_ext.replace('.', '_') # Remember the old filename, then give domain a # name with the extension old_domain_filename = self.domain.get_name() if not self.recursion: self.domain.set_name(old_domain_filename+filename_ext) # Change the domain starttime to the current time old_domain_starttime = self.domain.starttime self.domain.starttime = self.domain.time # Build a new data_structure. next_data_structure=\ Data_format_sww(self.domain, mode=self.mode,\ max_size = self.max_size,\ recursion = self.recursion+1) if not self.recursion: print ' file_size = %s'%file_size print ' saving file to %s'%next_data_structure.filename #set up the new data_structure self.domain.writer = next_data_structure #FIXME - could be cleaner to use domain.store_timestep etc. next_data_structure.store_connectivity() next_data_structure.store_timestep(names) fid.sync() fid.close() #restore the old starttime and filename self.domain.starttime = old_domain_starttime self.domain.set_name(old_domain_filename) else: self.recursion = False 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) if type(names) not in [types.ListType, types.TupleType]: names = [names] if 'stage' in names and 'xmomentum' in names and \ 'ymomentum' in names: # Get stage Q = domain.quantities['stage'] A,_ = Q.get_vertex_values(xy = False, precision = self.precision) z = fid.variables['elevation'] stage = choose(A-z[:] >= self.minimum_storable_height, (z[:], A)) # Get xmomentum Q = domain.quantities['xmomentum'] xmomentum, _ = Q.get_vertex_values(xy = False, precision = self.precision) # Get ymomentum Q = domain.quantities['ymomentum'] ymomentum, _ = Q.get_vertex_values(xy = False, precision = self.precision) # Write quantities to NetCDF self.writer.store_quantities(fid, time=self.domain.time, precision=self.precision, stage=stage, xmomentum=xmomentum, ymomentum=ymomentum) else: # This is producing a sww that is not standard. # Store time time[i] = self.domain.time 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': z = fid.variables['elevation'] A = choose(A-z[:] >= self.minimum_storable_height, (z[:], A)) 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 (Ole): But we need a UNIT test for that before # refactoring # Update extrema if requested domain = self.domain if domain.quantities_to_be_monitored is not None: for q, info in domain.quantities_to_be_monitored.items(): if info['min'] is not None: fid.variables[q + '.extrema'][0] = info['min'] fid.variables[q + '.min_location'][:] =\ info['min_location'] fid.variables[q + '.min_time'][0] = info['min_time'] if info['max'] is not None: fid.variables[q + '.extrema'][1] = info['max'] fid.variables[q + '.max_location'][:] =\ info['max_location'] fid.variables[q + '.max_time'][0] = info['max_time'] #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 DataFileNotOPenError, 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() #### NED is national exposure database (name changed to NEXIS) LAT_TITLE = 'LATITUDE' LONG_TITLE = 'LONGITUDE' X_TITLE = 'x' Y_TITLE = 'y' class Exposure_csv: def __init__(self,file_name, latitude_title=LAT_TITLE, longitude_title=LONG_TITLE, is_x_y_locations=None, x_title=X_TITLE, y_title=Y_TITLE, refine_polygon=None, title_check_list=None): """ This class is for handling the exposure csv file. It reads the file in and converts the lats and longs to a geospatial data object. Use the methods to read and write columns. The format of the csv files it reads is; The first row is a title row. comma's are the delimiters each column is a 'set' of data Feel free to use/expand it to read other csv files. It is not for adding and deleting rows Can geospatial handle string attributes? It's not made for them. Currently it can't load and save string att's. So just use geospatial to hold the x, y and georef? Bad, since different att's are in diferent structures. Not so bad, the info to write if the .csv file is saved is in attribute_dic The location info is in the geospatial attribute. """ self._file_name = file_name self._geospatial = None # # self._attribute_dic is a dictionary. #The keys are the column titles. #The values are lists of column data # self._title_index_dic is a dictionary. #The keys are the column titles. #The values are the index positions of file columns. self._attribute_dic, self._title_index_dic = \ csv2dict(self._file_name, title_check_list=title_check_list) try: #Have code here that handles caps or lower lats = self._attribute_dic[latitude_title] longs = self._attribute_dic[longitude_title] except KeyError: # maybe a warning.. #Let's see if this works.. if False != is_x_y_locations: is_x_y_locations = True pass else: self._geospatial = Geospatial_data(latitudes = lats, longitudes = longs) if is_x_y_locations is True: if self._geospatial is not None: pass #fixme throw an error try: xs = self._attribute_dic[x_title] ys = self._attribute_dic[y_title] points = [[float(i),float(j)] for i,j in map(None,xs,ys)] except KeyError: # maybe a warning.. msg = "Could not find location information." raise TitleValueError, msg else: self._geospatial = Geospatial_data(data_points=points) # create a list of points that are in the refining_polygon # described by a list of indexes representing the points def __cmp__(self, other): #print "self._attribute_dic",self._attribute_dic #print "other._attribute_dic",other._attribute_dic #print "self._title_index_dic", self._title_index_dic #print "other._title_index_dic", other._title_index_dic #check that a is an instance of this class if isinstance(self, type(other)): result = cmp(self._attribute_dic, other._attribute_dic) if result <>0: return result # The order of the columns is important. Therefore.. result = cmp(self._title_index_dic, other._title_index_dic) if result <>0: return result for self_ls, other_ls in map(None,self._attribute_dic, \ other._attribute_dic): result = cmp(self._attribute_dic[self_ls], other._attribute_dic[other_ls]) if result <>0: return result return 0 else: return 1 def get_column(self, column_name, use_refind_polygon=False): """ Given a column name return a list of the column values Note, the type of the values will be String! do this to change a list of strings to a list of floats time = [float(x) for x in time] Not implemented: if use_refind_polygon is True, only return values in the refined polygon """ if not self._attribute_dic.has_key(column_name): msg = 'Therer is no column called %s!' %column_name raise TitleValueError, msg return self._attribute_dic[column_name] def get_value(self, value_column_name, known_column_name, known_values, use_refind_polygon=False): """ Do linear interpolation on the known_colum, using the known_value, to return a value of the column_value_name. """ pass def get_location(self, use_refind_polygon=False): """ Return a geospatial object which describes the locations of the location file. Note, if there is not location info, this returns None. Not implemented: if use_refind_polygon is True, only return values in the refined polygon """ return self._geospatial def set_column(self, column_name, column_values, overwrite=False): """ Add a column to the 'end' (with the right most column being the end) of the csv file. Set overwrite to True if you want to overwrite a column. Note, in column_name white space is removed and case is not checked. Precondition The column_name and column_values cannot have comma's in it. """ value_row_count = \ len(self._attribute_dic[self._title_index_dic.keys()[0]]) if len(column_values) <> value_row_count: msg = 'The number of column values must equal the number of rows.' raise DataMissingValuesError, msg if self._attribute_dic.has_key(column_name): if not overwrite: msg = 'Column name %s already in use!' %column_name raise TitleValueError, msg else: # New title. Add it to the title index. self._title_index_dic[column_name] = len(self._title_index_dic) self._attribute_dic[column_name] = column_values #print "self._title_index_dic[column_name]",self._title_index_dic def save(self, file_name=None): """ Save the exposure csv file """ if file_name is None: file_name = self._file_name fd = open(file_name,'wb') writer = csv.writer(fd) #Write the title to a cvs file line = [None]* len(self._title_index_dic) for title in self._title_index_dic.iterkeys(): line[self._title_index_dic[title]]= title writer.writerow(line) # Write the values to a cvs file value_row_count = \ len(self._attribute_dic[self._title_index_dic.keys()[0]]) for row_i in range(value_row_count): line = [None]* len(self._title_index_dic) for title in self._title_index_dic.iterkeys(): line[self._title_index_dic[title]]= \ self._attribute_dic[title][row_i] writer.writerow(line) def csv2dict(file_name, title_check_list=None): """ Load in the csv as a dic, title as key and column info as value, . Also, create a dic, title as key and column index as value, to keep track of the column order. Two dictionaries are returned. WARNING: Vaules are returned as strings. do this to change a list of strings to a list of floats time = [float(x) for x in time] """ # attribute_dic = {} title_index_dic = {} titles_stripped = [] # list of titles reader = csv.reader(file(file_name)) # Read in and manipulate the title info titles = reader.next() for i,title in enumerate(titles): titles_stripped.append(title.strip()) title_index_dic[title.strip()] = i title_count = len(titles_stripped) #print "title_index_dic",title_index_dic if title_check_list is not None: for title_check in title_check_list: #msg = "Reading error. This row is not present ", title_check #assert title_index_dic.has_key(title_check), msg if not title_index_dic.has_key(title_check): #reader.close() msg = "Reading error. This row is not present ", \ title_check raise IOError, msg #create a dic of colum values, indexed by column title for line in reader: if len(line) <> title_count: raise IOError #FIXME make this nicer for i, value in enumerate(line): attribute_dic.setdefault(titles_stripped[i],[]).append(value) return attribute_dic, title_index_dic #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 anuga.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 DataTimeError, '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, easting_min=None, easting_max=None, northing_min=None, northing_max=None, use_cache=False, 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 """ kwargs = {'basename_out': basename_out, 'easting_min': easting_min, 'easting_max': easting_max, 'northing_min': northing_min, 'northing_max': northing_max, 'verbose': verbose} if use_cache is True: from caching import cache result = cache(_dem2pts, basename_in, kwargs, dependencies = [basename_in + '.dem'], verbose = verbose) else: result = apply(_dem2pts, [basename_in], kwargs) return result def _dem2pts(basename_in, basename_out=None, verbose=False, easting_min=None, easting_max=None, northing_min=None, northing_max=None): """Read Digitial Elevation model from the following NetCDF format (.dem) Internal function. See public function dem2pts for details. """ # FIXME: Can this be written feasibly using write_pts? import os from Scientific.IO.NetCDF import NetCDFFile from Numeric import Float, zeros, reshape, sum 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' # Assign default values xllcorner = xllcorner + 0.5*cellsize # Convert to gridline registration yllcorner = yllcorner + 0.5*cellsize if easting_min is None: easting_min = xllcorner if easting_max is None: easting_max = xllcorner + (ncols-1)*cellsize if northing_min is None: northing_min = yllcorner if northing_max is None: northing_max = yllcorner + (nrows-1)*cellsize # Compute offsets to update georeferencing easting_offset = xllcorner - easting_min northing_offset = yllcorner - northing_min # Georeferencing outfile.zone = zone outfile.xllcorner = easting_min # Easting of lower left corner outfile.yllcorner = northing_min # 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 dem_elevation_r = reshape(dem_elevation, (nrows, ncols)) totalnopoints = nrows*ncols # Calculating number of NODATA_values for each row in clipped region # FIXME: use array operations to do faster nn = 0 k = 0 i1_0 = 0 j1_0 = 0 thisj = 0 thisi = 0 for i in range(nrows): y = (nrows-i-1)*cellsize + yllcorner for j in range(ncols): x = j*cellsize + xllcorner if easting_min <= x <= easting_max and \ northing_min <= y <= northing_max: thisj = j thisi = i if dem_elevation_r[i,j] == NODATA_value: nn += 1 if k == 0: i1_0 = i j1_0 = j k += 1 index1 = j1_0 index2 = thisj # Dimension definitions nrows_in_bounding_box = int(round((northing_max-northing_min)/cellsize)) ncols_in_bounding_box = int(round((easting_max-easting_min)/cellsize)) clippednopoints = (thisi+1-i1_0)*(thisj+1-j1_0) nopoints = clippednopoints-nn clipped_dem_elev = dem_elevation_r[i1_0:thisi+1,j1_0:thisj+1] if verbose: print 'There are %d values in the elevation' %totalnopoints print 'There are %d values in the clipped elevation' %clippednopoints print 'There are %d NODATA_values in the clipped elevation' %nn outfile.createDimension('number_of_points', nopoints) 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'] lenv = index2-index1+1 # Store data global_index = 0 # for i in range(nrows): for i in range(i1_0,thisi+1,1): if verbose and i%((nrows+10)/10)==0: print 'Processing row %d of %d' %(i, nrows) lower_index = global_index v = dem_elevation_r[i,index1:index2+1] no_NODATA = sum(v == NODATA_value) if no_NODATA > 0: newcols = lenv - no_NODATA # ncols_in_bounding_box - no_NODATA else: newcols = lenv # ncols_in_bounding_box telev = zeros(newcols, Float) tpoints = zeros((newcols, 2), Float) local_index = 0 y = (nrows-i-1)*cellsize + yllcorner #for j in range(ncols): for j in range(j1_0,index2+1,1): x = j*cellsize + xllcorner if easting_min <= x <= easting_max and \ northing_min <= y <= northing_max and \ dem_elevation_r[i,j] <> NODATA_value: tpoints[local_index, :] = [x-easting_min,y-northing_min] telev[local_index] = dem_elevation_r[i, j] global_index += 1 local_index += 1 upper_index = global_index if upper_index == lower_index + newcols: points[lower_index:upper_index, :] = tpoints elevation[lower_index:upper_index] = telev assert global_index == nopoints, 'index not equal to number of points' infile.close() outfile.close() def _read_hecras_cross_sections(lines): """Return block of surface lines for each cross section Starts with SURFACE LINE, Ends with END CROSS-SECTION """ points = [] reading_surface = False for i, line in enumerate(lines): if len(line.strip()) == 0: #Ignore blanks continue if lines[i].strip().startswith('SURFACE LINE'): reading_surface = True continue if lines[i].strip().startswith('END') and reading_surface: yield points reading_surface = False points = [] if reading_surface: fields = line.strip().split(',') easting = float(fields[0]) northing = float(fields[1]) elevation = float(fields[2]) points.append([easting, northing, elevation]) def hecras_cross_sections2pts(basename_in, basename_out=None, verbose=False): """Read HEC-RAS Elevation datal from the following ASCII format (.sdf) Example: # RAS export file created on Mon 15Aug2005 11:42 # by HEC-RAS Version 3.1.1 BEGIN HEADER: UNITS: METRIC DTM TYPE: TIN DTM: v:\1\cit\perth_topo\river_tin STREAM LAYER: c:\local\hecras\21_02_03\up_canning_cent3d.shp CROSS-SECTION LAYER: c:\local\hecras\21_02_03\up_can_xs3d.shp MAP PROJECTION: UTM PROJECTION ZONE: 50 DATUM: AGD66 VERTICAL DATUM: NUMBER OF REACHES: 19 NUMBER OF CROSS-SECTIONS: 14206 END HEADER: Only the SURFACE LINE data of the following form will be utilised CROSS-SECTION: STREAM ID:Southern-Wungong REACH ID:Southern-Wungong STATION:19040.* CUT LINE: 405548.671603161 , 6438142.7594925 405734.536092045 , 6438326.10404912 405745.130459356 , 6438331.48627354 405813.89633823 , 6438368.6272789 SURFACE LINE: 405548.67, 6438142.76, 35.37 405552.24, 6438146.28, 35.41 405554.78, 6438148.78, 35.44 405555.80, 6438149.79, 35.44 405559.37, 6438153.31, 35.45 405560.88, 6438154.81, 35.44 405562.93, 6438156.83, 35.42 405566.50, 6438160.35, 35.38 405566.99, 6438160.83, 35.37 ... END CROSS-SECTION 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, zeros, reshape root = basename_in #Get ASCII file infile = open(root + '.sdf', 'r') #Open SDF file for read if verbose: print 'Reading DEM from %s' %(root + '.sdf') lines = infile.readlines() infile.close() if verbose: print 'Converting to pts format' i = 0 while lines[i].strip() == '' or lines[i].strip().startswith('#'): i += 1 assert lines[i].strip().upper() == 'BEGIN HEADER:' i += 1 assert lines[i].strip().upper().startswith('UNITS:') units = lines[i].strip().split()[1] i += 1 assert lines[i].strip().upper().startswith('DTM TYPE:') i += 1 assert lines[i].strip().upper().startswith('DTM:') i += 1 assert lines[i].strip().upper().startswith('STREAM') i += 1 assert lines[i].strip().upper().startswith('CROSS') i += 1 assert lines[i].strip().upper().startswith('MAP PROJECTION:') projection = lines[i].strip().split(':')[1] i += 1 assert lines[i].strip().upper().startswith('PROJECTION ZONE:') zone = int(lines[i].strip().split(':')[1]) i += 1 assert lines[i].strip().upper().startswith('DATUM:') datum = lines[i].strip().split(':')[1] i += 1 assert lines[i].strip().upper().startswith('VERTICAL DATUM:') i += 1 assert lines[i].strip().upper().startswith('NUMBER OF REACHES:') i += 1 assert lines[i].strip().upper().startswith('NUMBER OF CROSS-SECTIONS:') number_of_cross_sections = int(lines[i].strip().split(':')[1]) i += 1 #Now read all points points = [] elevation = [] for j, entries in enumerate(_read_hecras_cross_sections(lines[i:])): for k, entry in enumerate(entries): points.append(entry[:2]) elevation.append(entry[2]) msg = 'Actual #number_of_cross_sections == %d, Reported as %d'\ %(j+1, number_of_cross_sections) assert j+1 == number_of_cross_sections, msg #Get output file if basename_out == None: ptsname = root + '.pts' else: ptsname = basename_out + '.pts' geo_ref = Geo_reference(zone, 0, 0, datum, projection, units) geo = Geospatial_data(points, {"elevation":elevation}, verbose=verbose, geo_reference=geo_ref) geo.export_points_file(ptsname) def export_grid(basename_in, extra_name_out = None, quantities = None, # defaults to elevation timestep = None, reduction = None, cellsize = 10, NODATA_value = -9999, easting_min = None, easting_max = None, northing_min = None, northing_max = None, verbose = False, origin = None, datum = 'WGS84', format = 'ers'): """ Wrapper for sww2dem. - see sww2dem to find out what most of the parameters do. Quantities is a list of quantities. Each quantity will be calculated for each sww file. This returns the basenames of the files returned, which is made up of the dir and all of the file name, except the extension. This function returns the names of the files produced. It will also produce as many output files as there are input sww files. """ if quantities is None: quantities = ['elevation'] if type(quantities) is str: quantities = [quantities] # How many sww files are there? dir, base = os.path.split(basename_in) #print "basename_in",basename_in #print "base",base iterate_over = get_all_swwfiles(dir,base,verbose) if dir == "": dir = "." # Unix compatibility files_out = [] #print 'sww_file',iterate_over for sww_file in iterate_over: for quantity in quantities: if extra_name_out is None: basename_out = sww_file + '_' + quantity else: basename_out = sww_file + '_' + quantity + '_' \ + extra_name_out # print "basename_out", basename_out file_out = sww2dem(dir+sep+sww_file, dir+sep+basename_out, quantity, timestep, reduction, cellsize, NODATA_value, easting_min, easting_max, northing_min, northing_max, verbose, origin, datum, format) files_out.append(file_out) #print "basenames_out after",basenames_out return files_out def get_timeseries(production_dirs, output_dir, scenario_name, gauges_dir_name, plot_quantity, generate_fig = False, reportname = None, surface = False, time_min = None, time_max = None, title_on = False, verbose = True, nodes=None): """ nodes - number of processes used. warning - this function has no tests """ if reportname == None: report = False else: report = True if nodes is None: is_parallel = False else: is_parallel = True # Generate figures swwfiles = {} if is_parallel is True: for i in range(nodes): print 'Sending node %d of %d' %(i,nodes) swwfiles = {} if not reportname == None: reportname = report_name + '_%s' %(i) for label_id in production_dirs.keys(): if label_id == 'boundaries': swwfile = best_boundary_sww else: file_loc = output_dir + label_id + sep sww_extra = '_P%s_%s' %(i,nodes) swwfile = file_loc + scenario_name + sww_extra + '.sww' print 'swwfile',swwfile swwfiles[swwfile] = label_id texname, elev_output = sww2timeseries(swwfiles, gauges_dir_name, production_dirs, report = report, reportname = reportname, plot_quantity = plot_quantity, generate_fig = generate_fig, surface = surface, time_min = time_min, time_max = time_max, title_on = title_on, verbose = verbose) else: for label_id in production_dirs.keys(): if label_id == 'boundaries': print 'boundaries' file_loc = project.boundaries_in_dir swwfile = project.boundaries_dir_name3 + '.sww' # swwfile = boundary_dir_filename else: file_loc = output_dir + label_id + sep swwfile = file_loc + scenario_name + '.sww' swwfiles[swwfile] = label_id texname, elev_output = sww2timeseries(swwfiles, gauges_dir_name, production_dirs, report = report, reportname = reportname, plot_quantity = plot_quantity, generate_fig = generate_fig, surface = surface, time_min = time_min, time_max = time_max, title_on = title_on, verbose = verbose) def sww2dem(basename_in, basename_out = None, quantity = None, # defaults to elevation timestep = None, reduction = None, cellsize = 10, NODATA_value = -9999, easting_min = None, easting_max = None, northing_min = None, northing_max = None, verbose = False, origin = None, datum = 'WGS84', format = 'ers'): """Read SWW file and convert to Digitial Elevation model format (.asc or .ers) Example (ASC): ncols 3121 nrows 1800 xllcorner 722000 yllcorner 5893000 cellsize 25 NODATA_value -9999 138.3698 137.4194 136.5062 135.5558 .......... Also write accompanying file with same basename_in but extension .prj 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 The parameter quantity must be the name of an existing quantity or an expression involving existing quantities. The default is 'elevation'. Quantity is not a list of quantities. if timestep (an index) is given, output quantity at that timestep if reduction is given use that to reduce quantity over all timesteps. datum format can be either 'asc' or 'ers' """ import sys from Numeric import array, Float, concatenate, NewAxis, zeros, reshape, \ sometrue from Numeric import array2string from anuga.utilities.polygon import inside_polygon, outside_polygon, \ separate_points_by_polygon from anuga.abstract_2d_finite_volumes.util import \ apply_expression_to_dictionary msg = 'Format must be either asc or ers' assert format.lower() in ['asc', 'ers'], msg false_easting = 500000 false_northing = 10000000 if quantity is None: quantity = 'elevation' if reduction is None: reduction = max if basename_out is None: basename_out = basename_in + '_%s' %quantity if quantity_formula.has_key(quantity): quantity = quantity_formula[quantity] swwfile = basename_in + '.sww' demfile = basename_out + '.' + format # Note the use of a .ers extension is optional (write_ermapper_grid will # deal with either option #if verbose: bye= nsuadsfd[0] # uncomment to check catching verbose errors # Read sww file if verbose: print 'Reading from %s' %swwfile print 'Output directory is %s' %basename_out from Scientific.IO.NetCDF import NetCDFFile fid = NetCDFFile(swwfile) #Get extent and reference x = fid.variables['x'][:] y = fid.variables['y'][:] volumes = fid.variables['volumes'][:] if timestep is not None: times = fid.variables['time'][timestep] else: times = fid.variables['time'][:] number_of_timesteps = fid.dimensions['number_of_timesteps'] number_of_points = fid.dimensions['number_of_points'] if origin is None: # Get geo_reference # sww files don't have to have a geo_ref try: geo_reference = Geo_reference(NetCDFObject=fid) except AttributeError, e: geo_reference = Geo_reference() # Default georef object xllcorner = geo_reference.get_xllcorner() yllcorner = geo_reference.get_yllcorner() zone = geo_reference.get_zone() else: zone = origin[0] xllcorner = origin[1] yllcorner = origin[2] # FIXME: Refactor using code from Interpolation_function.statistics # (in interpolate.py) # Something like print swwstats(swwname) if verbose: print '------------------------------------------------' print 'Statistics of SWW file:' print ' Name: %s' %swwfile print ' Reference:' print ' Lower left corner: [%f, %f]'\ %(xllcorner, yllcorner) if timestep is not None: print ' Time: %f' %(times) else: print ' Start time: %f' %fid.starttime[0] print ' Extent:' print ' x [m] in [%f, %f], len(x) == %d'\ %(min(x.flat), max(x.flat), len(x.flat)) print ' y [m] in [%f, %f], len(y) == %d'\ %(min(y.flat), max(y.flat), len(y.flat)) if timestep is not None: print ' t [s] = %f, len(t) == %d' %(times, 1) else: print ' t [s] in [%f, %f], len(t) == %d'\ %(min(times), max(times), len(times)) print ' Quantities [SI units]:' for name in ['stage', 'xmomentum', 'ymomentum']: q = fid.variables[name][:].flat if timestep is not None: q = q[timestep*len(x):(timestep+1)*len(x)] if verbose: print ' %s in [%f, %f]' %(name, min(q), max(q)) for name in ['elevation']: q = fid.variables[name][:].flat if verbose: print ' %s in [%f, %f]' %(name, min(q), max(q)) # Get quantity and reduce if applicable if verbose: print 'Processing quantity %s' %quantity # Turn NetCDF objects into Numeric arrays quantity_dict = {} for name in fid.variables.keys(): quantity_dict[name] = fid.variables[name][:] # Convert quantity expression to quantities found in sww file q = apply_expression_to_dictionary(quantity, quantity_dict) if len(q.shape) == 2: #q has a time component and needs to be reduced along #the temporal dimension if verbose: print 'Reducing quantity %s' %quantity q_reduced = zeros( number_of_points, Float ) if timestep is not None: for k in range(number_of_points): q_reduced[k] = q[timestep,k] else: for k in range(number_of_points): q_reduced[k] = reduction( q[:,k] ) q = q_reduced #Post condition: Now q has dimension: number_of_points assert len(q.shape) == 1 assert q.shape[0] == number_of_points if verbose: print 'Processed values for %s are in [%f, %f]' %(quantity, min(q), max(q)) #Create grid and update xll/yll corner and x,y #Relative extent if easting_min is None: xmin = min(x) else: xmin = easting_min - xllcorner if easting_max is None: xmax = max(x) else: xmax = easting_max - xllcorner if northing_min is None: ymin = min(y) else: ymin = northing_min - yllcorner if northing_max is None: ymax = max(y) else: ymax = northing_max - yllcorner if verbose: print 'Creating grid' ncols = int((xmax-xmin)/cellsize)+1 nrows = int((ymax-ymin)/cellsize)+1 #New absolute reference and coordinates newxllcorner = xmin+xllcorner newyllcorner = ymin+yllcorner x = x+xllcorner-newxllcorner y = y+yllcorner-newyllcorner vertex_points = concatenate ((x[:, NewAxis] ,y[:, NewAxis]), axis = 1) assert len(vertex_points.shape) == 2 grid_points = zeros ( (ncols*nrows, 2), Float ) for i in xrange(nrows): if format.lower() == 'asc': yg = i*cellsize else: #this will flip the order of the y values for ers yg = (nrows-i)*cellsize for j in xrange(ncols): xg = j*cellsize k = i*ncols + j grid_points[k,0] = xg grid_points[k,1] = yg #Interpolate from anuga.fit_interpolate.interpolate import Interpolate # Remove loners from vertex_points, volumes here vertex_points, volumes = remove_lone_verts(vertex_points, volumes) #export_mesh_file('monkey.tsh',{'vertices':vertex_points, 'triangles':volumes}) #import sys; sys.exit() interp = Interpolate(vertex_points, volumes, verbose = verbose) #Interpolate using quantity values if verbose: print 'Interpolating' grid_values = interp.interpolate(q, grid_points).flat if verbose: print 'Interpolated values are in [%f, %f]' %(min(grid_values), max(grid_values)) #Assign NODATA_value to all points outside bounding polygon (from interpolation mesh) P = interp.mesh.get_boundary_polygon() outside_indices = outside_polygon(grid_points, P, closed=True) for i in outside_indices: grid_values[i] = NODATA_value if format.lower() == 'ers': # setup ERS header information grid_values = reshape(grid_values,(nrows, ncols)) header = {} header['datum'] = '"' + datum + '"' # FIXME The use of hardwired UTM and zone number needs to be made optional # FIXME Also need an automatic test for coordinate type (i.e. EN or LL) header['projection'] = '"UTM-' + str(zone) + '"' header['coordinatetype'] = 'EN' if header['coordinatetype'] == 'LL': header['longitude'] = str(newxllcorner) header['latitude'] = str(newyllcorner) elif header['coordinatetype'] == 'EN': header['eastings'] = str(newxllcorner) header['northings'] = str(newyllcorner) header['nullcellvalue'] = str(NODATA_value) header['xdimension'] = str(cellsize) header['ydimension'] = str(cellsize) header['value'] = '"' + quantity + '"' #header['celltype'] = 'IEEE8ByteReal' #FIXME: Breaks unit test #Write if verbose: print 'Writing %s' %demfile import ermapper_grids ermapper_grids.write_ermapper_grid(demfile, grid_values, header) fid.close() else: #Write to Ascii format #Write prj file prjfile = basename_out + '.prj' if verbose: print 'Writing %s' %prjfile prjid = open(prjfile, 'w') prjid.write('Projection %s\n' %'UTM') prjid.write('Zone %d\n' %zone) prjid.write('Datum %s\n' %datum) prjid.write('Zunits NO\n') prjid.write('Units METERS\n') prjid.write('Spheroid %s\n' %datum) prjid.write('Xshift %d\n' %false_easting) prjid.write('Yshift %d\n' %false_northing) prjid.write('Parameters\n') prjid.close() if verbose: print 'Writing %s' %demfile ascid = open(demfile, 'w') ascid.write('ncols %d\n' %ncols) ascid.write('nrows %d\n' %nrows) ascid.write('xllcorner %d\n' %newxllcorner) ascid.write('yllcorner %d\n' %newyllcorner) ascid.write('cellsize %f\n' %cellsize) ascid.write('NODATA_value %d\n' %NODATA_value) #Get bounding polygon from mesh #P = interp.mesh.get_boundary_polygon() #inside_indices = inside_polygon(grid_points, P) for i in range(nrows): if verbose and i%((nrows+10)/10)==0: print 'Doing row %d of %d' %(i, nrows) base_index = (nrows-i-1)*ncols slice = grid_values[base_index:base_index+ncols] s = array2string(slice, max_line_width=sys.maxint) ascid.write(s[1:-1] + '\n') #print #for j in range(ncols): # index = base_index+j# # print grid_values[index], # ascid.write('%f ' %grid_values[index]) #ascid.write('\n') #Close ascid.close() fid.close() return basename_out #Backwards compatibility def sww2asc(basename_in, basename_out = None, quantity = None, timestep = None, reduction = None, cellsize = 10, verbose = False, origin = None): print 'sww2asc will soon be obsoleted - please use sww2dem' sww2dem(basename_in, basename_out = basename_out, quantity = quantity, timestep = timestep, reduction = reduction, cellsize = cellsize, verbose = verbose, origin = origin, datum = 'WGS84', format = 'asc') def sww2ers(basename_in, basename_out = None, quantity = None, timestep = None, reduction = None, cellsize = 10, verbose = False, origin = None, datum = 'WGS84'): print 'sww2ers will soon be obsoleted - please use sww2dem' sww2dem(basename_in, basename_out = basename_out, quantity = quantity, timestep = timestep, reduction = reduction, cellsize = cellsize, verbose = verbose, origin = origin, datum = datum, format = 'ers') ################################# END COMPATIBILITY ############## def sww2pts(basename_in, basename_out=None, data_points=None, quantity=None, timestep=None, reduction=None, NODATA_value=-9999, verbose=False, origin=None): #datum = 'WGS84') """Read SWW file and convert to interpolated values at selected points The parameter quantity' must be the name of an existing quantity or an expression involving existing quantities. The default is 'elevation'. if timestep (an index) is given, output quantity at that timestep if reduction is given use that to reduce quantity over all timesteps. data_points (Nx2 array) give locations of points where quantity is to be computed. """ import sys from Numeric import array, Float, concatenate, NewAxis, zeros, reshape, sometrue from Numeric import array2string from anuga.utilities.polygon import inside_polygon, outside_polygon, separate_points_by_polygon from anuga.abstract_2d_finite_volumes.util import apply_expression_to_dictionary from anuga.geospatial_data.geospatial_data import Geospatial_data if quantity is None: quantity = 'elevation' if reduction is None: reduction = max if basename_out is None: basename_out = basename_in + '_%s' %quantity swwfile = basename_in + '.sww' ptsfile = basename_out + '.pts' # Read sww file if verbose: print 'Reading from %s' %swwfile from Scientific.IO.NetCDF import NetCDFFile fid = NetCDFFile(swwfile) # Get extent and reference x = fid.variables['x'][:] y = fid.variables['y'][:] volumes = fid.variables['volumes'][:] number_of_timesteps = fid.dimensions['number_of_timesteps'] number_of_points = fid.dimensions['number_of_points'] if origin is None: # Get geo_reference # sww files don't have to have a geo_ref try: geo_reference = Geo_reference(NetCDFObject=fid) except AttributeError, e: geo_reference = Geo_reference() #Default georef object xllcorner = geo_reference.get_xllcorner() yllcorner = geo_reference.get_yllcorner() zone = geo_reference.get_zone() else: zone = origin[0] xllcorner = origin[1] yllcorner = origin[2] # FIXME: Refactor using code from file_function.statistics # Something like print swwstats(swwname) if verbose: x = fid.variables['x'][:] y = fid.variables['y'][:] times = fid.variables['time'][:] print '------------------------------------------------' print 'Statistics of SWW file:' print ' Name: %s' %swwfile print ' Reference:' print ' Lower left corner: [%f, %f]'\ %(xllcorner, yllcorner) print ' Start time: %f' %fid.starttime[0] print ' Extent:' print ' x [m] in [%f, %f], len(x) == %d'\ %(min(x.flat), max(x.flat), len(x.flat)) print ' y [m] in [%f, %f], len(y) == %d'\ %(min(y.flat), max(y.flat), len(y.flat)) print ' t [s] in [%f, %f], len(t) == %d'\ %(min(times), max(times), len(times)) print ' Quantities [SI units]:' for name in ['stage', 'xmomentum', 'ymomentum', 'elevation']: q = fid.variables[name][:].flat print ' %s in [%f, %f]' %(name, min(q), max(q)) # Get quantity and reduce if applicable if verbose: print 'Processing quantity %s' %quantity # Turn NetCDF objects into Numeric arrays quantity_dict = {} for name in fid.variables.keys(): quantity_dict[name] = fid.variables[name][:] # Convert quantity expression to quantities found in sww file q = apply_expression_to_dictionary(quantity, quantity_dict) if len(q.shape) == 2: # q has a time component and needs to be reduced along # the temporal dimension if verbose: print 'Reducing quantity %s' %quantity q_reduced = zeros( number_of_points, Float ) for k in range(number_of_points): q_reduced[k] = reduction( q[:,k] ) q = q_reduced # Post condition: Now q has dimension: number_of_points assert len(q.shape) == 1 assert q.shape[0] == number_of_points if verbose: print 'Processed values for %s are in [%f, %f]' %(quantity, min(q), max(q)) # Create grid and update xll/yll corner and x,y vertex_points = concatenate ((x[:, NewAxis] ,y[:, NewAxis]), axis = 1) assert len(vertex_points.shape) == 2 # Interpolate from anuga.fit_interpolate.interpolate import Interpolate interp = Interpolate(vertex_points, volumes, verbose = verbose) # Interpolate using quantity values if verbose: print 'Interpolating' interpolated_values = interp.interpolate(q, data_points).flat if verbose: print 'Interpolated values are in [%f, %f]' %(min(interpolated_values), max(interpolated_values)) # Assign NODATA_value to all points outside bounding polygon (from interpolation mesh) P = interp.mesh.get_boundary_polygon() outside_indices = outside_polygon(data_points, P, closed=True) for i in outside_indices: interpolated_values[i] = NODATA_value # Store results G = Geospatial_data(data_points=data_points, attributes=interpolated_values) G.export_points_file(ptsfile, absolute = True) fid.close() def convert_dem_from_ascii2netcdf(basename_in, basename_out = None, use_cache = False, 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 """ kwargs = {'basename_out': basename_out, 'verbose': verbose} if use_cache is True: from caching import cache result = cache(_convert_dem_from_ascii2netcdf, basename_in, kwargs, dependencies = [basename_in + '.asc', basename_in + '.prj'], verbose = verbose) else: result = apply(_convert_dem_from_ascii2netcdf, [basename_in], kwargs) return result def _convert_dem_from_ascii2netcdf(basename_in, basename_out = None, verbose = False): """Read Digitial Elevation model from the following ASCII format (.asc) Internal function. See public function convert_dem_from_ascii2netcdf for details. """ 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 n = len(lines[6:]) for i, line in enumerate(lines[6:]): fields = line.split() if verbose and i%((n+10)/10)==0: 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, fail_on_NaN = True, NaN_filler = 0, elevation = None, inverted_bathymetry = True ): #FIXME: Bathymetry should be obtained #from MOST somehow. #Alternatively from elsewhere #or, as a last resort, #specified here. #The value of -100 will work #for the Wollongong tsunami #scenario but is very hacky """Convert MOST and 'Ferret' NetCDF format for wave propagation to sww format native to abstract_2d_finite_volumes. 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 from Numeric import allclose, around precision = Float msg = 'Must use latitudes and longitudes for minlat, maxlon etc' if minlat != None: assert -90 < minlat < 90 , msg if maxlat != None: assert -90 < maxlat < 90 , msg if minlat != None: assert maxlat > minlat if minlon != None: assert -180 < minlon < 180 , msg if maxlon != None: assert -180 < maxlon < 180 , msg if minlon != None: assert maxlon > minlon #Get NetCDF data if verbose: print 'Reading files %s_*.nc' %basename_in #print "basename_in + '_ha.nc'",basename_in + '_ha.nc' 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) file_e = NetCDFFile(basename_in + '_e.nc', 'r') #Elevation (z) (m) if basename_out is None: swwname = basename_in + '.sww' else: swwname = basename_out + '.sww' # Get dimensions of file_h for dimension in file_h.dimensions.keys(): if dimension[:3] == 'LON': dim_h_longitude = dimension if dimension[:3] == 'LAT': dim_h_latitude = dimension if dimension[:4] == 'TIME': dim_h_time = dimension # print 'long:', dim_h_longitude # print 'lats:', dim_h_latitude # print 'times:', dim_h_time times = file_h.variables[dim_h_time] latitudes = file_h.variables[dim_h_latitude] longitudes = file_h.variables[dim_h_longitude] # get dimensions for file_e for dimension in file_e.dimensions.keys(): if dimension[:3] == 'LON': dim_e_longitude = dimension if dimension[:3] == 'LAT': dim_e_latitude = dimension # get dimensions for file_u for dimension in file_u.dimensions.keys(): if dimension[:3] == 'LON': dim_u_longitude = dimension if dimension[:3] == 'LAT': dim_u_latitude = dimension if dimension[:4] == 'TIME': dim_u_time = dimension # get dimensions for file_v for dimension in file_v.dimensions.keys(): if dimension[:3] == 'LON': dim_v_longitude = dimension if dimension[:3] == 'LAT': dim_v_latitude = dimension if dimension[:4] == 'TIME': dim_v_time = dimension # Precision used by most for lat/lon is 4 or 5 decimals e_lat = around(file_e.variables[dim_e_latitude][:], 5) e_lon = around(file_e.variables[dim_e_longitude][:], 5) # Check that files are compatible assert allclose(latitudes, file_u.variables[dim_u_latitude]) assert allclose(latitudes, file_v.variables[dim_v_latitude]) assert allclose(latitudes, e_lat) assert allclose(longitudes, file_u.variables[dim_u_longitude]) assert allclose(longitudes, file_v.variables[dim_v_longitude]) assert allclose(longitudes, e_lon) if mint is None: jmin = 0 mint = times[0] else: jmin = searchsorted(times, mint) if maxt is None: jmax = len(times) maxt = times[-1] else: jmax = searchsorted(times, maxt) #print "latitudes[:]",latitudes[:] #print "longitudes[:]",longitudes [:] kmin, kmax, lmin, lmax = _get_min_max_indexes(latitudes[:], longitudes[:], minlat, maxlat, minlon, maxlon) times = times[jmin:jmax] latitudes = latitudes[kmin:kmax] longitudes = longitudes[lmin:lmax] #print "latitudes[:]",latitudes[:] #print "longitudes[:]",longitudes [:] if verbose: print 'cropping' zname = 'ELEVATION' amplitudes = file_h.variables['HA'][jmin:jmax, kmin:kmax, lmin:lmax] uspeed = file_u.variables['UA'][jmin:jmax, kmin:kmax, lmin:lmax] #Lon vspeed = file_v.variables['VA'][jmin:jmax, kmin:kmax, lmin:lmax] #Lat elevations = file_e.variables[zname][kmin:kmax, lmin:lmax] # if latitudes2[0]==latitudes[0] and latitudes2[-1]==latitudes[-1]: # elevations = file_e.variables['ELEVATION'][kmin:kmax, lmin:lmax] # elif latitudes2[0]==latitudes[-1] and latitudes2[-1]==latitudes[0]: # from Numeric import asarray # elevations=elevations.tolist() # elevations.reverse() # elevations=asarray(elevations) # else: # from Numeric import asarray # elevations=elevations.tolist() # elevations.reverse() # elevations=asarray(elevations) # 'print hmmm' #Get missing values nan_ha = file_h.variables['HA'].missing_value[0] nan_ua = file_u.variables['UA'].missing_value[0] nan_va = file_v.variables['VA'].missing_value[0] if hasattr(file_e.variables[zname],'missing_value'): nan_e = file_e.variables[zname].missing_value[0] else: nan_e = None #Cleanup from Numeric import sometrue missing = (amplitudes == nan_ha) if sometrue (missing): if fail_on_NaN: msg = 'NetCDFFile %s contains missing values'\ %(basename_in+'_ha.nc') raise DataMissingValuesError, msg else: amplitudes = amplitudes*(missing==0) + missing*NaN_filler missing = (uspeed == nan_ua) if sometrue (missing): if fail_on_NaN: msg = 'NetCDFFile %s contains missing values'\ %(basename_in+'_ua.nc') raise DataMissingValuesError, msg else: uspeed = uspeed*(missing==0) + missing*NaN_filler missing = (vspeed == nan_va) if sometrue (missing): if fail_on_NaN: msg = 'NetCDFFile %s contains missing values'\ %(basename_in+'_va.nc') raise DataMissingValuesError, msg else: vspeed = vspeed*(missing==0) + missing*NaN_filler missing = (elevations == nan_e) if sometrue (missing): if fail_on_NaN: msg = 'NetCDFFile %s contains missing values'\ %(basename_in+'_e.nc') raise DataMissingValuesError, msg else: elevations = elevations*(missing==0) + missing*NaN_filler ####### 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 if verbose: print '------------------------------------------------' print 'Statistics:' print ' Extent (lat/lon):' print ' lat in [%f, %f], len(lat) == %d'\ %(min(latitudes.flat), max(latitudes.flat), len(latitudes.flat)) print ' lon in [%f, %f], len(lon) == %d'\ %(min(longitudes.flat), max(longitudes.flat), len(longitudes.flat)) print ' t in [%f, %f], len(t) == %d'\ %(min(times.flat), max(times.flat), len(times.flat)) q = amplitudes.flat name = 'Amplitudes (ha) [cm]' print ' %s in [%f, %f]' %(name, min(q), max(q)) q = uspeed.flat name = 'Speeds (ua) [cm/s]' print ' %s in [%f, %f]' %(name, min(q), max(q)) q = vspeed.flat name = 'Speeds (va) [cm/s]' print ' %s in [%f, %f]' %(name, min(q), max(q)) q = elevations.flat name = 'Elevations (e) [m]' print ' %s in [%f, %f]' %(name, min(q), max(q)) # 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 file_h.close() file_u.close() file_v.close() file_e.close() # NetCDF file definition outfile = NetCDFFile(swwname, 'w') description = 'Converted from Ferret files: %s, %s, %s, %s'\ %(basename_in + '_ha.nc', basename_in + '_ua.nc', basename_in + '_va.nc', basename_in + '_e.nc') # Create new file starttime = times[0] sww = Write_sww() sww.store_header(outfile, times, number_of_volumes, number_of_points, description=description, verbose=verbose) # Store from anuga.coordinate_transforms.redfearn import redfearn x = zeros(number_of_points, Float) #Easting y = zeros(number_of_points, Float) #Northing if verbose: print 'Making triangular grid' # Check zone boundaries refzone, _, _ = redfearn(latitudes[0],longitudes[0]) vertices = {} i = 0 for k, lat in enumerate(latitudes): #Y direction for l, lon in enumerate(longitudes): #X direction 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): #X direction for k in range(number_of_latitudes-1): #Y direction 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 is None: origin = Geo_reference(refzone,min(x),min(y)) geo_ref = write_NetCDF_georeference(origin, outfile) if elevation is not None: z = elevation else: if inverted_bathymetry: z = -1*elevations else: z = elevations #FIXME: z should be obtained from MOST and passed in here from Numeric import resize z = resize(z,outfile.variables['z'][:].shape) outfile.variables['x'][:] = x - geo_ref.get_xllcorner() outfile.variables['y'][:] = y - geo_ref.get_yllcorner() outfile.variables['z'][:] = z #FIXME HACK for bacwards compat. outfile.variables['elevation'][:] = z outfile.variables['volumes'][:] = volumes.astype(Int32) #For Opteron 64 #Time stepping stage = outfile.variables['stage'] xmomentum = outfile.variables['xmomentum'] ymomentum = outfile.variables['ymomentum'] if verbose: print 'Converting quantities' n = len(times) for j in range(n): if verbose and j%((n+10)/10)==0: print ' Doing %d of %d' %(j, n) i = 0 for k in range(number_of_latitudes): #Y direction for l in range(number_of_longitudes): #X direction w = zscale*amplitudes[j,k,l]/100 + mean_stage stage[j,i] = w h = w - z[i] xmomentum[j,i] = uspeed[j,k,l]/100*h ymomentum[j,i] = vspeed[j,k,l]/100*h i += 1 #outfile.close() #FIXME: Refactor using code from file_function.statistics #Something like print swwstats(swwname) if verbose: x = outfile.variables['x'][:] y = outfile.variables['y'][:] print '------------------------------------------------' print 'Statistics of output file:' print ' Name: %s' %swwname print ' Reference:' print ' Lower left corner: [%f, %f]'\ %(geo_ref.get_xllcorner(), geo_ref.get_yllcorner()) print ' Start time: %f' %starttime print ' Min time: %f' %mint print ' Max time: %f' %maxt print ' Extent:' print ' x [m] in [%f, %f], len(x) == %d'\ %(min(x.flat), max(x.flat), len(x.flat)) print ' y [m] in [%f, %f], len(y) == %d'\ %(min(y.flat), max(y.flat), len(y.flat)) print ' t [s] in [%f, %f], len(t) == %d'\ %(min(times), max(times), len(times)) print ' Quantities [SI units]:' for name in ['stage', 'xmomentum', 'ymomentum', 'elevation']: q = outfile.variables[name][:].flat print ' %s in [%f, %f]' %(name, min(q), max(q)) outfile.close() def timefile2netcdf(filename, quantity_names=None, time_as_seconds=False): """Template for converting typical text files with time series to NetCDF tms file. The file format is assumed to be either two fields separated by a comma: time [DD/MM/YY hh:mm:ss], value0 value1 value2 ... E.g 31/08/04 00:00:00, 1.328223 0 0 31/08/04 00:15:00, 1.292912 0 0 or time (seconds), value0 value1 value2 ... 0.0, 1.328223 0 0 0.1, 1.292912 0 0 will provide a time dependent function f(t) with three attributes filename is assumed to be the rootname with extenisons .txt and .sww """ import time, calendar from Numeric import array from anuga.config import time_format from anuga.utilities.numerical_tools import ensure_numeric fid = open(filename + '.txt') line = fid.readline() fid.close() fields = line.split(',') msg = 'File %s must have the format date, value0 value1 value2 ...' assert len(fields) == 2, msg if not time_as_seconds: try: starttime = calendar.timegm(time.strptime(fields[0], time_format)) except ValueError: msg = 'First field in file %s must be' %filename msg += ' date-time with format %s.\n' %time_format msg += 'I got %s instead.' %fields[0] raise DataTimeError, msg else: try: starttime = float(fields[0]) except Error: msg = "Bad time format" raise DataTimeError, msg #Split values values = [] for value in fields[1].split(): values.append(float(value)) q = ensure_numeric(values) msg = 'ERROR: File must contain at least one independent value' assert len(q.shape) == 1, msg #Read times proper from Numeric import zeros, Float, alltrue from anuga.config import time_format import time, calendar fid = open(filename + '.txt') lines = fid.readlines() fid.close() N = len(lines) d = len(q) T = zeros(N, Float) #Time Q = zeros((N, d), Float) #Values for i, line in enumerate(lines): fields = line.split(',') if not time_as_seconds: realtime = calendar.timegm(time.strptime(fields[0], time_format)) else: realtime = float(fields[0]) T[i] = realtime - starttime for j, value in enumerate(fields[1].split()): Q[i, j] = float(value) msg = 'File %s must list time as a monotonuosly ' %filename msg += 'increasing sequence' assert alltrue( T[1:] - T[:-1] > 0 ), msg #Create NetCDF file from Scientific.IO.NetCDF import NetCDFFile fid = NetCDFFile(filename + '.tms', 'w') fid.institution = 'Geoscience Australia' fid.description = 'Time series' #Reference point #Start time in seconds since the epoch (midnight 1/1/1970) #FIXME: Use Georef fid.starttime = starttime # dimension definitions #fid.createDimension('number_of_volumes', self.number_of_volumes) #fid.createDimension('number_of_vertices', 3) fid.createDimension('number_of_timesteps', len(T)) fid.createVariable('time', Float, ('number_of_timesteps',)) fid.variables['time'][:] = T for i in range(Q.shape[1]): try: name = quantity_names[i] except: name = 'Attribute%d'%i fid.createVariable(name, Float, ('number_of_timesteps',)) fid.variables[name][:] = Q[:,i] fid.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)] def sww2domain(filename,boundary=None,t=None,\ fail_if_NaN=True,NaN_filler=0\ ,verbose = False,very_verbose = False): """ Usage: domain = sww2domain('file.sww',t=time (default = last time in file)) Boundary is not recommended if domain.smooth is not selected, as it uses unique coordinates, but not unique boundaries. This means that the boundary file will not be compatable with the coordinates, and will give a different final boundary, or crash. """ NaN=9.969209968386869e+036 #initialise NaN. from Scientific.IO.NetCDF import NetCDFFile from shallow_water import Domain from Numeric import asarray, transpose, resize if verbose: print 'Reading from ', filename fid = NetCDFFile(filename, 'r') #Open existing file for read time = fid.variables['time'] #Timesteps if t is None: t = time[-1] time_interp = get_time_interp(time,t) # Get the variables as Numeric arrays x = fid.variables['x'][:] #x-coordinates of vertices y = fid.variables['y'][:] #y-coordinates of vertices elevation = fid.variables['elevation'] #Elevation stage = fid.variables['stage'] #Water level xmomentum = fid.variables['xmomentum'] #Momentum in the x-direction ymomentum = fid.variables['ymomentum'] #Momentum in the y-direction starttime = fid.starttime[0] volumes = fid.variables['volumes'][:] #Connectivity coordinates=transpose(asarray([x.tolist(),y.tolist()])) conserved_quantities = [] interpolated_quantities = {} other_quantities = [] # get geo_reference #sww files don't have to have a geo_ref try: geo_reference = Geo_reference(NetCDFObject=fid) except: #AttributeError, e: geo_reference = None if verbose: print ' getting quantities' for quantity in fid.variables.keys(): dimensions = fid.variables[quantity].dimensions if 'number_of_timesteps' in dimensions: conserved_quantities.append(quantity) interpolated_quantities[quantity]=\ interpolated_quantity(fid.variables[quantity][:],time_interp) else: other_quantities.append(quantity) other_quantities.remove('x') other_quantities.remove('y') other_quantities.remove('z') other_quantities.remove('volumes') try: other_quantities.remove('stage_range') other_quantities.remove('xmomentum_range') other_quantities.remove('ymomentum_range') other_quantities.remove('elevation_range') except: pass conserved_quantities.remove('time') if verbose: print ' building domain' # From domain.Domain: # domain = Domain(coordinates, volumes,\ # conserved_quantities = conserved_quantities,\ # other_quantities = other_quantities,zone=zone,\ # xllcorner=xllcorner, yllcorner=yllcorner) # From shallow_water.Domain: coordinates=coordinates.tolist() volumes=volumes.tolist() #FIXME:should this be in mesh?(peter row) if fid.smoothing == 'Yes': unique = False else: unique = True if unique: coordinates,volumes,boundary=weed(coordinates,volumes,boundary) try: domain = Domain(coordinates, volumes, boundary) except AssertionError, e: fid.close() msg = 'Domain could not be created: %s. Perhaps use "fail_if_NaN=False and NaN_filler = ..."' %e raise DataDomainError, msg if not boundary is None: domain.boundary = boundary domain.geo_reference = geo_reference domain.starttime=float(starttime)+float(t) domain.time=0.0 for quantity in other_quantities: try: NaN = fid.variables[quantity].missing_value except: pass #quantity has no missing_value number X = fid.variables[quantity][:] if very_verbose: print ' ',quantity print ' NaN =',NaN print ' max(X)' print ' ',max(X) print ' max(X)==NaN' print ' ',max(X)==NaN print '' if (max(X)==NaN) or (min(X)==NaN): if fail_if_NaN: msg = 'quantity "%s" contains no_data entry'%quantity raise DataMissingValuesError, msg else: data = (X<>NaN) X = (X*data)+(data==0)*NaN_filler if unique: X = resize(X,(len(X)/3,3)) domain.set_quantity(quantity,X) # for quantity in conserved_quantities: try: NaN = fid.variables[quantity].missing_value except: pass #quantity has no missing_value number X = interpolated_quantities[quantity] if very_verbose: print ' ',quantity print ' NaN =',NaN print ' max(X)' print ' ',max(X) print ' max(X)==NaN' print ' ',max(X)==NaN print '' if (max(X)==NaN) or (min(X)==NaN): if fail_if_NaN: msg = 'quantity "%s" contains no_data entry'%quantity raise DataMissingValuesError, msg else: data = (X<>NaN) X = (X*data)+(data==0)*NaN_filler if unique: X = resize(X,(X.shape[0]/3,3)) domain.set_quantity(quantity,X) fid.close() return domain def interpolated_quantity(saved_quantity,time_interp): #given an index and ratio, interpolate quantity with respect to time. index,ratio = time_interp Q = saved_quantity if ratio > 0: q = (1-ratio)*Q[index]+ ratio*Q[index+1] else: q = Q[index] #Return vector of interpolated values return q def get_time_interp(time,t=None): #Finds the ratio and index for time interpolation. #It is borrowed from previous abstract_2d_finite_volumes code. if t is None: t=time[-1] index = -1 ratio = 0. else: T = time tau = t index=0 msg = 'Time interval derived from file %s [%s:%s]'\ %('FIXMEfilename', T[0], T[-1]) msg += ' does not match model time: %s' %tau if tau < time[0]: raise DataTimeError, msg if tau > time[-1]: raise DataTimeError, msg while tau > time[index]: index += 1 while tau < time[index]: index -= 1 if tau == time[index]: #Protect against case where tau == time[-1] (last time) # - also works in general when tau == time[i] ratio = 0 else: #t is now between index and index+1 ratio = (tau - time[index])/(time[index+1] - time[index]) return (index,ratio) def weed(coordinates,volumes,boundary = None): if type(coordinates)==ArrayType: coordinates = coordinates.tolist() if type(volumes)==ArrayType: volumes = volumes.tolist() unique = False point_dict = {} same_point = {} for i in range(len(coordinates)): point = tuple(coordinates[i]) if point_dict.has_key(point): unique = True same_point[i]=point #to change all point i references to point j else: point_dict[point]=i same_point[i]=point coordinates = [] i = 0 for point in point_dict.keys(): point = tuple(point) coordinates.append(list(point)) point_dict[point]=i i+=1 for volume in volumes: for i in range(len(volume)): index = volume[i] if index>-1: volume[i]=point_dict[same_point[index]] new_boundary = {} if not boundary is None: for segment in boundary.keys(): point0 = point_dict[same_point[segment[0]]] point1 = point_dict[same_point[segment[1]]] label = boundary[segment] #FIXME should the bounday attributes be concaterated #('exterior, pond') or replaced ('pond')(peter row) if new_boundary.has_key((point0,point1)): new_boundary[(point0,point1)]=new_boundary[(point0,point1)]#\ #+','+label elif new_boundary.has_key((point1,point0)): new_boundary[(point1,point0)]=new_boundary[(point1,point0)]#\ #+','+label else: new_boundary[(point0,point1)]=label boundary = new_boundary return coordinates,volumes,boundary def decimate_dem(basename_in, stencil, cellsize_new, 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 .......... Decimate data to cellsize_new using stencil and write to NetCDF dem format. """ import os from Scientific.IO.NetCDF import NetCDFFile from Numeric import Float, zeros, sum, reshape, equal root = basename_in inname = root + '.dem' #Open existing netcdf file to read infile = NetCDFFile(inname, 'r') if verbose: print 'Reading DEM from %s' %inname #Read metadata ncols = infile.ncols[0] nrows = infile.nrows[0] xllcorner = infile.xllcorner[0] yllcorner = infile.yllcorner[0] cellsize = infile.cellsize[0] NODATA_value = infile.NODATA_value[0] zone = infile.zone[0] false_easting = infile.false_easting[0] false_northing = infile.false_northing[0] projection = infile.projection datum = infile.datum units = infile.units dem_elevation = infile.variables['elevation'] #Get output file name if basename_out == None: outname = root + '_' + repr(cellsize_new) + '.dem' else: outname = basename_out + '.dem' if verbose: print 'Write decimated NetCDF file to %s' %outname #Determine some dimensions for decimated grid (nrows_stencil, ncols_stencil) = stencil.shape x_offset = ncols_stencil / 2 y_offset = nrows_stencil / 2 cellsize_ratio = int(cellsize_new / cellsize) ncols_new = 1 + (ncols - ncols_stencil) / cellsize_ratio nrows_new = 1 + (nrows - nrows_stencil) / cellsize_ratio #Open netcdf file for output outfile = NetCDFFile(outname, 'w') #Create new file outfile.institution = 'Geoscience Australia' outfile.description = 'NetCDF DEM format for compact and portable storage ' +\ 'of spatial point data' #Georeferencing outfile.zone = zone outfile.projection = projection outfile.datum = datum outfile.units = units outfile.cellsize = cellsize_new outfile.NODATA_value = NODATA_value outfile.false_easting = false_easting outfile.false_northing = false_northing outfile.xllcorner = xllcorner + (x_offset * cellsize) outfile.yllcorner = yllcorner + (y_offset * cellsize) outfile.ncols = ncols_new outfile.nrows = nrows_new # dimension definition outfile.createDimension('number_of_points', nrows_new*ncols_new) # variable definition outfile.createVariable('elevation', Float, ('number_of_points',)) # Get handle to the variable elevation = outfile.variables['elevation'] dem_elevation_r = reshape(dem_elevation, (nrows, ncols)) #Store data global_index = 0 for i in range(nrows_new): if verbose: print 'Processing row %d of %d' %(i, nrows_new) lower_index = global_index telev = zeros(ncols_new, Float) local_index = 0 trow = i * cellsize_ratio for j in range(ncols_new): tcol = j * cellsize_ratio tmp = dem_elevation_r[trow:trow+nrows_stencil, tcol:tcol+ncols_stencil] #if dem contains 1 or more NODATA_values set value in #decimated dem to NODATA_value, else compute decimated #value using stencil if sum(sum(equal(tmp, NODATA_value))) > 0: telev[local_index] = NODATA_value else: telev[local_index] = sum(sum(tmp * stencil)) global_index += 1 local_index += 1 upper_index = global_index elevation[lower_index:upper_index] = telev assert global_index == nrows_new*ncols_new, 'index not equal to number of points' infile.close() outfile.close() def tsh2sww(filename, verbose=False): #test_tsh2sww """ to check if a tsh/msh file 'looks' good. """ if verbose == True:print 'Creating domain from', filename domain = pmesh_to_domain_instance(filename, Domain) if verbose == True:print "Number of triangles = ", len(domain) domain.smooth = True domain.format = 'sww' #Native netcdf visualisation format file_path, filename = path.split(filename) filename, ext = path.splitext(filename) domain.set_name(filename) domain.reduction = mean if verbose == True:print "file_path",file_path if file_path == "":file_path = "." domain.set_datadir(file_path) if verbose == True: print "Output written to " + domain.get_datadir() + sep + \ domain.get_name() + "." + domain.format sww = get_dataobject(domain) sww.store_connectivity() sww.store_timestep('stage') def asc_csiro2sww(bath_dir, elevation_dir, ucur_dir, vcur_dir, sww_file, minlat = None, maxlat = None, minlon = None, maxlon = None, zscale=1, mean_stage = 0, fail_on_NaN = True, elevation_NaN_filler = 0, bath_prefix='ba', elevation_prefix='el', verbose=False): """ Produce an sww boundary file, from esri ascii data from CSIRO. Also convert latitude and longitude to UTM. All coordinates are assumed to be given in the GDA94 datum. assume: All files are in esri ascii format 4 types of information bathymetry elevation u velocity v velocity Assumptions The metadata of all the files is the same Each type is in a seperate directory One bath file with extention .000 The time period is less than 24hrs and uniform. """ from Scientific.IO.NetCDF import NetCDFFile from anuga.coordinate_transforms.redfearn import redfearn precision = Float # So if we want to change the precision its done here # go in to the bath dir and load the only file, bath_files = os.listdir(bath_dir) bath_file = bath_files[0] bath_dir_file = bath_dir + os.sep + bath_file bath_metadata,bath_grid = _read_asc(bath_dir_file) #Use the date.time of the bath file as a basis for #the start time for other files base_start = bath_file[-12:] #go into the elevation dir and load the 000 file elevation_dir_file = elevation_dir + os.sep + elevation_prefix \ + base_start elevation_files = os.listdir(elevation_dir) ucur_files = os.listdir(ucur_dir) vcur_files = os.listdir(vcur_dir) elevation_files.sort() # the first elevation file should be the # file with the same base name as the bath data assert elevation_files[0] == 'el' + base_start number_of_latitudes = bath_grid.shape[0] number_of_longitudes = bath_grid.shape[1] number_of_volumes = (number_of_latitudes-1)*(number_of_longitudes-1)*2 longitudes = [bath_metadata['xllcorner']+x*bath_metadata['cellsize'] \ for x in range(number_of_longitudes)] latitudes = [bath_metadata['yllcorner']+y*bath_metadata['cellsize'] \ for y in range(number_of_latitudes)] # reverse order of lat, so the fist lat represents the first grid row latitudes.reverse() kmin, kmax, lmin, lmax = _get_min_max_indexes(latitudes[:],longitudes[:], minlat=minlat, maxlat=maxlat, minlon=minlon, maxlon=maxlon) bath_grid = bath_grid[kmin:kmax,lmin:lmax] latitudes = latitudes[kmin:kmax] longitudes = longitudes[lmin:lmax] number_of_latitudes = len(latitudes) number_of_longitudes = len(longitudes) number_of_times = len(os.listdir(elevation_dir)) number_of_points = number_of_latitudes*number_of_longitudes number_of_volumes = (number_of_latitudes-1)*(number_of_longitudes-1)*2 #Work out the times if len(elevation_files) > 1: # Assume: The time period is less than 24hrs. time_period = (int(elevation_files[1][-3:]) - \ int(elevation_files[0][-3:]))*60*60 times = [x*time_period for x in range(len(elevation_files))] else: times = [0.0] if verbose: print '------------------------------------------------' print 'Statistics:' print ' Extent (lat/lon):' print ' lat in [%f, %f], len(lat) == %d'\ %(min(latitudes), max(latitudes), len(latitudes)) print ' lon in [%f, %f], len(lon) == %d'\ %(min(longitudes), max(longitudes), len(longitudes)) print ' t in [%f, %f], len(t) == %d'\ %(min(times), max(times), len(times)) ######### WRITE THE SWW FILE ############# # NetCDF file definition outfile = NetCDFFile(sww_file, 'w') #Create new file outfile.institution = 'Geoscience Australia' outfile.description = 'Converted from XXX' #For sww compatibility outfile.smoothing = 'Yes' outfile.order = 1 #Start time in seconds since the epoch (midnight 1/1/1970) outfile.starttime = 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', number_of_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 anuga.coordinate_transforms.redfearn import redfearn x = zeros(number_of_points, Float) #Easting y = zeros(number_of_points, Float) #Northing if verbose: print 'Making triangular grid' #Get zone of 1st point. refzone, _, _ = redfearn(latitudes[0],longitudes[0]) vertices = {} i = 0 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): #X direction for k in range(number_of_latitudes-1): #Y direction v1 = vertices[l,k+1] v2 = vertices[l,k] v3 = vertices[l+1,k+1] v4 = vertices[l+1,k] #Note, this is different to the ferrit2sww code #since the order of the lats is reversed. volumes.append([v1,v3,v2]) #Upper element volumes.append([v4,v2,v3]) #Lower element volumes = array(volumes) geo_ref = Geo_reference(refzone,min(x),min(y)) geo_ref.write_NetCDF(outfile) # This will put the geo ref in the middle #geo_ref = Geo_reference(refzone,(max(x)+min(x))/2.0,(max(x)+min(y))/2.) if verbose: print '------------------------------------------------' print 'More Statistics:' print ' Extent (/lon):' print ' x in [%f, %f], len(lat) == %d'\ %(min(x), max(x), len(x)) print ' y in [%f, %f], len(lon) == %d'\ %(min(y), max(y), len(y)) print 'geo_ref: ',geo_ref z = resize(bath_grid,outfile.variables['z'][:].shape) outfile.variables['x'][:] = x - geo_ref.get_xllcorner() outfile.variables['y'][:] = y - geo_ref.get_yllcorner() outfile.variables['z'][:] = z outfile.variables['elevation'][:] = z #FIXME HACK outfile.variables['volumes'][:] = volumes.astype(Int32) #On Opteron 64 stage = outfile.variables['stage'] xmomentum = outfile.variables['xmomentum'] ymomentum = outfile.variables['ymomentum'] outfile.variables['time'][:] = times #Store time relative if verbose: print 'Converting quantities' n = number_of_times for j in range(number_of_times): # load in files elevation_meta, elevation_grid = \ _read_asc(elevation_dir + os.sep + elevation_files[j]) _, u_momentum_grid = _read_asc(ucur_dir + os.sep + ucur_files[j]) _, v_momentum_grid = _read_asc(vcur_dir + os.sep + vcur_files[j]) #cut matrix to desired size elevation_grid = elevation_grid[kmin:kmax,lmin:lmax] u_momentum_grid = u_momentum_grid[kmin:kmax,lmin:lmax] v_momentum_grid = v_momentum_grid[kmin:kmax,lmin:lmax] # handle missing values missing = (elevation_grid == elevation_meta['NODATA_value']) if sometrue (missing): if fail_on_NaN: msg = 'File %s contains missing values'\ %(elevation_files[j]) raise DataMissingValuesError, msg else: elevation_grid = elevation_grid*(missing==0) + \ missing*elevation_NaN_filler if verbose and j%((n+10)/10)==0: print ' Doing %d of %d' %(j, n) i = 0 for k in range(number_of_latitudes): #Y direction for l in range(number_of_longitudes): #X direction w = zscale*elevation_grid[k,l] + mean_stage stage[j,i] = w h = w - z[i] xmomentum[j,i] = u_momentum_grid[k,l]*h ymomentum[j,i] = v_momentum_grid[k,l]*h i += 1 outfile.close() def _get_min_max_indexes(latitudes_ref,longitudes_ref, minlat=None, maxlat=None, minlon=None, maxlon=None): """ return max, min indexes (for slicing) of the lat and long arrays to cover the area specified with min/max lat/long Think of the latitudes and longitudes describing a 2d surface. The area returned is, if possible, just big enough to cover the inputed max/min area. (This will not be possible if the max/min area has a section outside of the latitudes/longitudes area.) asset longitudes are sorted, long - from low to high (west to east, eg 148 - 151) assert latitudes are sorted, ascending or decending """ latitudes = latitudes_ref[:] longitudes = longitudes_ref[:] latitudes = ensure_numeric(latitudes) longitudes = ensure_numeric(longitudes) assert allclose(sort(longitudes), longitudes) lat_ascending = True if not allclose(sort(latitudes), latitudes): lat_ascending = False # reverse order of lat, so it's in ascending order latitudes = latitudes[::-1] assert allclose(sort(latitudes), latitudes) #print "latitudes in funct", latitudes largest_lat_index = len(latitudes)-1 #Cut out a smaller extent. if minlat == None: lat_min_index = 0 else: lat_min_index = searchsorted(latitudes, minlat)-1 if lat_min_index <0: lat_min_index = 0 if maxlat == None: lat_max_index = largest_lat_index #len(latitudes) else: lat_max_index = searchsorted(latitudes, maxlat) if lat_max_index > largest_lat_index: lat_max_index = largest_lat_index if minlon == None: lon_min_index = 0 else: lon_min_index = searchsorted(longitudes, minlon)-1 if lon_min_index <0: lon_min_index = 0 if maxlon == None: lon_max_index = len(longitudes) else: lon_max_index = searchsorted(longitudes, maxlon) # Reversing the indexes, if the lat array is decending if lat_ascending is False: lat_min_index, lat_max_index = largest_lat_index - lat_max_index , \ largest_lat_index - lat_min_index lat_max_index = lat_max_index + 1 # taking into account how slicing works lon_max_index = lon_max_index + 1 # taking into account how slicing works return lat_min_index, lat_max_index, lon_min_index, lon_max_index def _read_asc(filename, verbose=False): """Read esri file 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 .......... """ datafile = open(filename) if verbose: print 'Reading DEM from %s' %(filename) lines = datafile.readlines() datafile.close() if verbose: print 'Got', len(lines), ' lines' ncols = int(lines.pop(0).split()[1].strip()) nrows = int(lines.pop(0).split()[1].strip()) xllcorner = float(lines.pop(0).split()[1].strip()) yllcorner = float(lines.pop(0).split()[1].strip()) cellsize = float(lines.pop(0).split()[1].strip()) NODATA_value = float(lines.pop(0).split()[1].strip()) assert len(lines) == nrows #Store data grid = [] n = len(lines) for i, line in enumerate(lines): cells = line.split() assert len(cells) == ncols grid.append(array([float(x) for x in cells])) grid = array(grid) return {'xllcorner':xllcorner, 'yllcorner':yllcorner, 'cellsize':cellsize, 'NODATA_value':NODATA_value}, grid #### URS 2 SWW ### lon_name = 'LON' lat_name = 'LAT' time_name = 'TIME' precision = Float # So if we want to change the precision its done here class Write_nc: """ Write an nc file. Note, this should be checked to meet cdc netcdf conventions for gridded data. http://www.cdc.noaa.gov/cdc/conventions/cdc_netcdf_standard.shtml """ def __init__(self, quantity_name, file_name, time_step_count, time_step, lon, lat): """ time_step_count is the number of time steps. time_step is the time step size pre-condition: quantity_name must be 'HA' 'UA'or 'VA'. """ self.quantity_name = quantity_name quantity_units = {'HA':'CENTIMETERS', 'UA':'CENTIMETERS/SECOND', 'VA':'CENTIMETERS/SECOND'} multiplier_dic = {'HA':100.0, # To convert from m to cm 'UA':100.0, # m/s to cm/sec 'VA':-100.0} # MUX files have positve x in the # Southern direction. This corrects for it, when writing nc files. self.quantity_multiplier = multiplier_dic[self.quantity_name] #self.file_name = file_name self.time_step_count = time_step_count self.time_step = time_step # NetCDF file definition self.outfile = NetCDFFile(file_name, 'w') outfile = self.outfile #Create new file nc_lon_lat_header(outfile, lon, lat) # TIME outfile.createDimension(time_name, None) outfile.createVariable(time_name, precision, (time_name,)) #QUANTITY outfile.createVariable(self.quantity_name, precision, (time_name, lat_name, lon_name)) outfile.variables[self.quantity_name].missing_value=-1.e+034 outfile.variables[self.quantity_name].units= \ quantity_units[self.quantity_name] outfile.variables[lon_name][:]= ensure_numeric(lon) outfile.variables[lat_name][:]= ensure_numeric(lat) #Assume no one will be wanting to read this, while we are writing #outfile.close() def store_timestep(self, quantity_slice): """ Write a time slice of quantity info quantity_slice is the data to be stored at this time step """ outfile = self.outfile # Get the variables time = outfile.variables[time_name] quantity = outfile.variables[self.quantity_name] i = len(time) #Store time time[i] = i*self.time_step #self.domain.time quantity[i,:] = quantity_slice* self.quantity_multiplier def close(self): self.outfile.close() def urs2sww(basename_in='o', basename_out=None, verbose=False, remove_nc_files=True, minlat=None, maxlat=None, minlon= None, maxlon=None, mint=None, maxt=None, mean_stage=0, origin = None, zscale=1, fail_on_NaN=True, NaN_filler=0, elevation=None): """ Convert URS C binary format for wave propagation to sww format native to abstract_2d_finite_volumes. Specify only basename_in and read files of the form basefilename_velocity-z-mux, basefilename_velocity-e-mux and basefilename_waveheight-n-mux 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. The latitude and longitude information is for a grid. 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 assumed to be in decimal degrees. NOTE: minlon is the most east boundary. origin is a 3-tuple with geo referenced UTM coordinates (zone, easting, northing) It will be the origin of the sww file. This shouldn't be used, since all of anuga should be able to handle an arbitary origin. URS C binary format has data orgainised as TIME, LONGITUDE, LATITUDE which means that latitude is the fastest varying dimension (row major order, so to speak) In URS C binary the latitudes and longitudes are in assending order. """ if basename_out == None: basename_out = basename_in files_out = urs2nc(basename_in, basename_out) ferret2sww(basename_out, minlat=minlat, maxlat=maxlat, minlon=minlon, maxlon=maxlon, mint=mint, maxt=maxt, mean_stage=mean_stage, origin=origin, zscale=zscale, fail_on_NaN=fail_on_NaN, NaN_filler=NaN_filler, inverted_bathymetry=True, verbose=verbose) #print "files_out",files_out if remove_nc_files: for file_out in files_out: os.remove(file_out) def urs2nc(basename_in = 'o', basename_out = 'urs'): """ Convert the 3 urs files to 4 nc files. The name of the urs file names must be; [basename_in]_velocity-z-mux [basename_in]_velocity-e-mux [basename_in]_waveheight-n-mux """ files_in = [basename_in + WAVEHEIGHT_MUX_LABEL, basename_in + EAST_VELOCITY_LABEL, basename_in + NORTH_VELOCITY_LABEL] files_out = [basename_out+'_ha.nc', basename_out+'_ua.nc', basename_out+'_va.nc'] quantities = ['HA','UA','VA'] for file_name in files_in: if os.access(file_name, os.F_OK) == 0 : msg = 'File %s does not exist or is not accessible' %file_name raise IOError, msg hashed_elevation = None for file_in, file_out, quantity in map(None, files_in, files_out, quantities): lonlatdep, lon, lat, depth = _binary_c2nc(file_in, file_out, quantity) #print "lonlatdep", lonlatdep if hashed_elevation == None: elevation_file = basename_out+'_e.nc' write_elevation_nc(elevation_file, lon, lat, depth) hashed_elevation = myhash(lonlatdep) else: msg = "The elevation information in the mux files is inconsistent" assert hashed_elevation == myhash(lonlatdep), msg files_out.append(elevation_file) return files_out def _binary_c2nc(file_in, file_out, quantity): """ Reads in a quantity urs file and writes a quantity nc file. additionally, returns the depth and lat, long info, so it can be written to a file. """ columns = 3 # long, lat , depth mux_file = open(file_in, 'rb') # Number of points/stations (points_num,)= unpack('i',mux_file.read(4)) # nt, int - Number of time steps (time_step_count,)= unpack('i',mux_file.read(4)) #dt, float - time step, seconds (time_step,) = unpack('f', mux_file.read(4)) msg = "Bad data in the mux file." if points_num < 0: mux_file.close() raise ANUGAError, msg if time_step_count < 0: mux_file.close() raise ANUGAError, msg if time_step < 0: mux_file.close() raise ANUGAError, msg lonlatdep = p_array.array('f') lonlatdep.read(mux_file, columns * points_num) lonlatdep = array(lonlatdep, typecode=Float) lonlatdep = reshape(lonlatdep, (points_num, columns)) lon, lat, depth = lon_lat2grid(lonlatdep) lon_sorted = list(lon) lon_sorted.sort() if not lon == lon_sorted: msg = "Longitudes in mux file are not in ascending order" raise IOError, msg lat_sorted = list(lat) lat_sorted.sort() if not lat == lat_sorted: msg = "Latitudes in mux file are not in ascending order" nc_file = Write_nc(quantity, file_out, time_step_count, time_step, lon, lat) for i in range(time_step_count): #Read in a time slice from mux file hz_p_array = p_array.array('f') hz_p_array.read(mux_file, points_num) hz_p = array(hz_p_array, typecode=Float) hz_p = reshape(hz_p, (len(lon), len(lat))) hz_p = transpose(hz_p) #mux has lat varying fastest, nc has long v.f. #write time slice to nc file nc_file.store_timestep(hz_p) mux_file.close() nc_file.close() return lonlatdep, lon, lat, depth def write_elevation_nc(file_out, lon, lat, depth_vector): """ Write an nc elevation file. """ # NetCDF file definition outfile = NetCDFFile(file_out, 'w') #Create new file nc_lon_lat_header(outfile, lon, lat) # ELEVATION zname = 'ELEVATION' outfile.createVariable(zname, precision, (lat_name, lon_name)) outfile.variables[zname].units='CENTIMETERS' outfile.variables[zname].missing_value=-1.e+034 outfile.variables[lon_name][:]= ensure_numeric(lon) outfile.variables[lat_name][:]= ensure_numeric(lat) depth = reshape(depth_vector, ( len(lat), len(lon))) outfile.variables[zname][:]= depth outfile.close() def nc_lon_lat_header(outfile, lon, lat): """ outfile is the netcdf file handle. lon - a list/array of the longitudes lat - a list/array of the latitudes """ outfile.institution = 'Geoscience Australia' outfile.description = 'Converted from URS binary C' # Longitude outfile.createDimension(lon_name, len(lon)) outfile.createVariable(lon_name, precision, (lon_name,)) outfile.variables[lon_name].point_spacing='uneven' outfile.variables[lon_name].units='degrees_east' outfile.variables[lon_name].assignValue(lon) # Latitude outfile.createDimension(lat_name, len(lat)) outfile.createVariable(lat_name, precision, (lat_name,)) outfile.variables[lat_name].point_spacing='uneven' outfile.variables[lat_name].units='degrees_north' outfile.variables[lat_name].assignValue(lat) def lon_lat2grid(long_lat_dep): """ given a list of points that are assumed to be an a grid, return the long's and lat's of the grid. long_lat_dep is an array where each row is a position. The first column is longitudes. The second column is latitudes. The latitude is the fastest varying dimension - in mux files """ LONG = 0 LAT = 1 QUANTITY = 2 long_lat_dep = ensure_numeric(long_lat_dep, Float) num_points = long_lat_dep.shape[0] this_rows_long = long_lat_dep[0,LONG] # Count the length of unique latitudes i = 0 while long_lat_dep[i,LONG] == this_rows_long and i < num_points: i += 1 # determine the lats and longsfrom the grid lat = long_lat_dep[:i, LAT] long = long_lat_dep[::i, LONG] lenlong = len(long) lenlat = len(lat) #print 'len lat', lat, len(lat) #print 'len long', long, len(long) msg = 'Input data is not gridded' assert num_points % lenlat == 0, msg assert num_points % lenlong == 0, msg # Test that data is gridded for i in range(lenlong): msg = 'Data is not gridded. It must be for this operation' first = i*lenlat last = first + lenlat assert allclose(long_lat_dep[first:last,LAT], lat), msg assert allclose(long_lat_dep[first:last,LONG], long[i]), msg # print 'range long', min(long), max(long) # print 'range lat', min(lat), max(lat) # print 'ref long', min(long_lat_dep[:,0]), max(long_lat_dep[:,0]) # print 'ref lat', min(long_lat_dep[:,1]), max(long_lat_dep[:,1]) msg = 'Out of range latitudes/longitudes' for l in lat:assert -90 < l < 90 , msg for l in long:assert -180 < l < 180 , msg # Changing quantity from lat being the fastest varying dimension to # long being the fastest varying dimension # FIXME - make this faster/do this a better way # use numeric transpose, after reshaping the quantity vector # quantity = zeros(len(long_lat_dep), Float) quantity = zeros(num_points, Float) # print 'num',num_points for lat_i, _ in enumerate(lat): for long_i, _ in enumerate(long): q_index = lat_i*lenlong+long_i lld_index = long_i*lenlat+lat_i # print 'lat_i', lat_i, 'long_i',long_i, 'q_index', q_index, 'lld_index', lld_index temp = long_lat_dep[lld_index, QUANTITY] quantity[q_index] = temp return long, lat, quantity #### END URS 2 SWW ### #### URS UNGRIDDED 2 SWW ### ### PRODUCING THE POINTS NEEDED FILE ### LL_LAT = -50.0 LL_LONG = 80.0 GRID_SPACING = 1.0/60.0 LAT_AMOUNT = 4800 LONG_AMOUNT = 3600 def URS_points_needed_to_file(file_name, boundary_polygon, zone, ll_lat=LL_LAT, ll_long=LL_LONG, grid_spacing=GRID_SPACING, lat_amount=LAT_AMOUNT, long_amount=LONG_AMOUNT, export_csv=False, use_cache=False, verbose=False): """ file_name - name of the urs file produced for David. boundary_polygon - a list of points that describes a polygon. The last point is assumed ot join the first point. This is in UTM (lat long would be better though) ll_lat - lower left latitude, in decimal degrees ll-long - lower left longitude, in decimal degrees grid_spacing - in deciamal degrees Don't add the file extension. It will be added. """ geo = URS_points_needed(boundary_polygon, zone, ll_lat, ll_long, grid_spacing, lat_amount, long_amount,use_cache, verbose) if not file_name[-4:] == ".urs": file_name += ".urs" geo.export_points_file(file_name) if export_csv: if file_name[-4:] == ".urs": file_name = file_name[:-4] + ".csv" geo.export_points_file(file_name) def URS_points_needed(boundary_polygon, zone, ll_lat=LL_LAT, ll_long=LL_LONG, grid_spacing=GRID_SPACING, lat_amount=LAT_AMOUNT, long_amount=LONG_AMOUNT, use_cache=False, verbose=False): args = (boundary_polygon, zone) kwargs = {'ll_lat': ll_lat, 'll_long': ll_long, 'grid_spacing': grid_spacing, 'lat_amount': lat_amount, 'long_amount': long_amount} if use_cache is True: try: from anuga.caching import cache except: msg = 'Caching was requested, but caching module'+\ 'could not be imported' raise msg geo = cache(_URS_points_needed, args, kwargs, verbose=verbose, compression=False) else: #I was getting 'got multiple values for keyword argument' errors #geo = apply(_URS_points_needed, args, kwargs) geo = _URS_points_needed(boundary_polygon, zone, ll_lat, ll_long, grid_spacing, lat_amount, long_amount) return geo def _URS_points_needed(boundary_polygon, zone, ll_lat=LL_LAT, ll_long=LL_LONG, grid_spacing=GRID_SPACING, lat_amount=LAT_AMOUNT, long_amount=LONG_AMOUNT): """ boundary_polygon - a list of points that describes a polygon. The last point is assumed ot join the first point. This is in UTM (lat long would be better though) ll_lat - lower left latitude, in decimal degrees ll-long - lower left longitude, in decimal degrees grid_spacing - in deciamal degrees """ from sets import ImmutableSet msg = "grid_spacing can not be zero" assert not grid_spacing ==0, msg a = boundary_polygon # List of segments. Each segment is two points. segs = [i and [a[i-1], a[i]] or [a[len(a)-1], a[0]] for i in range(len(a))] # convert the segs to Lat's and longs. # Don't assume the zone of the segments is the same as the lower left # corner of the lat long data!! They can easily be in different zones lat_long_set = ImmutableSet() for seg in segs: points_lat_long = points_needed(seg, ll_lat, ll_long, grid_spacing, lat_amount, long_amount, zone) lat_long_set |= ImmutableSet(points_lat_long) #print "lat_long_set",lat_long_set geo = Geospatial_data(data_points=list(lat_long_set), points_are_lats_longs=True) return geo def points_needed(seg, ll_lat, ll_long, grid_spacing, lat_amount, long_amount, zone): """ return a list of the points, in lats and longs that are needed to interpolate any point on the segment. """ from math import sqrt #print "zone",zone geo_reference = Geo_reference(zone=zone) #print "seg",seg geo = Geospatial_data(seg,geo_reference=geo_reference) seg_lat_long = geo.get_data_points(as_lat_long=True) #print "seg_lat_long", seg_lat_long # 1.415 = 2^0.5, rounded up.... sqrt_2_rounded_up = 1.415 buffer = sqrt_2_rounded_up * grid_spacing max_lat = max(seg_lat_long[0][0], seg_lat_long[1][0]) + buffer max_long = max(seg_lat_long[0][1], seg_lat_long[1][1]) + buffer min_lat = min(seg_lat_long[0][0], seg_lat_long[1][0]) - buffer min_long = min(seg_lat_long[0][1], seg_lat_long[1][1]) - buffer first_row = (min_long - ll_long)/grid_spacing # To round up first_row_long = int(round(first_row + 0.5)) #print "first_row", first_row_long last_row = (max_long - ll_long)/grid_spacing # round down last_row_long = int(round(last_row)) #print "last_row",last_row _long first_row = (min_lat - ll_lat)/grid_spacing # To round up first_row_lat = int(round(first_row + 0.5)) #print "first_row", first_row_lat last_row = (max_lat - ll_lat)/grid_spacing # round down last_row_lat = int(round(last_row)) #print "last_row",last_row_lat # to work out the max distance - # 111120 - horizontal distance between 1 deg latitude. #max_distance = sqrt_2_rounded_up * 111120 * grid_spacing max_distance = 157147.4112 * grid_spacing #print "max_distance", max_distance #2619.12 m for 1 minute points_lat_long = [] # Create a list of the lat long points to include. for index_lat in range(first_row_lat, last_row_lat + 1): for index_long in range(first_row_long, last_row_long + 1): lat = ll_lat + index_lat*grid_spacing long = ll_long + index_long*grid_spacing #filter here to keep good points if keep_point(lat, long, seg, max_distance): points_lat_long.append((lat, long)) #must be hashable #print "points_lat_long", points_lat_long # Now that we have these points, lets throw ones out that are too far away return points_lat_long def keep_point(lat, long, seg, max_distance): """ seg is two points, UTM """ from math import sqrt _ , x0, y0 = redfearn(lat, long) x1 = seg[0][0] y1 = seg[0][1] x2 = seg[1][0] y2 = seg[1][1] x2_1 = x2-x1 y2_1 = y2-y1 d = abs((x2_1)*(y1-y0)-(x1-x0)*(y2_1))/sqrt((x2_1)*(x2_1)+(y2_1)*(y2_1)) if d <= max_distance: return True else: return False #### CONVERTING UNGRIDDED URS DATA TO AN SWW FILE #### WAVEHEIGHT_MUX_LABEL = '_waveheight-z-mux' EAST_VELOCITY_LABEL = '_velocity-e-mux' NORTH_VELOCITY_LABEL = '_velocity-n-mux' def urs_ungridded2sww(basename_in='o', basename_out=None, verbose=False, mint=None, maxt=None, mean_stage=0, origin=None, hole_points_UTM=None, zscale=1): """ Convert URS C binary format for wave propagation to sww format native to abstract_2d_finite_volumes. Specify only basename_in and read files of the form basefilename_velocity-z-mux, basefilename_velocity-e-mux and basefilename_waveheight-n-mux 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. The latitude and longitude information is assumed ungridded grid. min's and max's: If omitted - full extend is used. To include a value min ans max may equal it. Lat and lon are assumed to be in decimal degrees. origin is a 3-tuple with geo referenced UTM coordinates (zone, easting, northing) It will be the origin of the sww file. This shouldn't be used, since all of anuga should be able to handle an arbitary origin. The mux point info is NOT relative to this origin. URS C binary format has data orgainised as TIME, LONGITUDE, LATITUDE which means that latitude is the fastest varying dimension (row major order, so to speak) In URS C binary the latitudes and longitudes are in assending order. Note, interpolations of the resulting sww file will be different from results of urs2sww. This is due to the interpolation function used, and the different grid structure between urs2sww and this function. Interpolating data that has an underlying gridded source can easily end up with different values, depending on the underlying mesh. consider these 4 points 50 -50 0 0 The grid can be - |\| A - or; - |/| B - If a point is just below the center of the midpoint, it will have a +ve value in grid A and a -ve value in grid B. """ from anuga.pmesh.mesh import Mesh, NoTrianglesError files_in = [basename_in + WAVEHEIGHT_MUX_LABEL, basename_in + EAST_VELOCITY_LABEL, basename_in + NORTH_VELOCITY_LABEL] quantities = ['HA','UA','VA'] # instanciate urs_points of the three mux files. mux = {} for quantity, file in map(None, quantities, files_in): mux[quantity] = Urs_points(file) # Could check that the depth is the same. (hashing) # handle to a mux file to do depth stuff a_mux = mux[quantities[0]] # Convert to utm lat = a_mux.lonlatdep[:,1] long = a_mux.lonlatdep[:,0] points_utm, zone = convert_from_latlon_to_utm( \ latitudes=lat, longitudes=long) #print "points_utm", points_utm #print "zone", zone elevation = a_mux.lonlatdep[:,2] * -1 # # grid ( create a mesh from the selected points) # This mesh has a problem. Triangles are streched over ungridded areas. # If these areas could be described as holes in pmesh, that would be great # I can't just get the user to selection a point in the middle. # A boundary is needed around these points. # But if the zone of points is obvious enough auto-segment should do # a good boundary. mesh = Mesh() mesh.add_vertices(points_utm) mesh.auto_segment(smooth_indents=True, expand_pinch=True) # To try and avoid alpha shape 'hugging' too much mesh.auto_segment( mesh.shape.get_alpha()*1.1 ) if hole_points_UTM is not None: point = ensure_absolute(hole_points_UTM) mesh.add_hole(point[0], point[1]) try: mesh.generate_mesh(minimum_triangle_angle=0.0, verbose=False) except NoTrianglesError: # This is a bit of a hack, going in and changing the # data structure. mesh.holes = [] mesh.generate_mesh(minimum_triangle_angle=0.0, verbose=False) mesh_dic = mesh.Mesh2MeshList() #mesh.export_mesh_file(basename_in + '_168.tsh') #import sys; sys.exit() # These are the times of the mux file mux_times = [] for i in range(a_mux.time_step_count): mux_times.append(a_mux.time_step * i) mux_times_start_i, mux_times_fin_i = mux2sww_time(mux_times, mint, maxt) times = mux_times[mux_times_start_i:mux_times_fin_i] if mux_times_start_i == mux_times_fin_i: # Close the mux files for quantity, file in map(None, quantities, files_in): mux[quantity].close() msg="Due to mint and maxt there's no time info in the boundary SWW." raise Exception, msg # If this raise is removed there is currently no downstream errors points_utm=ensure_numeric(points_utm) #print "mesh_dic['generatedpointlist']", mesh_dic['generatedpointlist'] #print "points_utm", points_utm assert ensure_numeric(mesh_dic['generatedpointlist']) == \ ensure_numeric(points_utm) volumes = mesh_dic['generatedtrianglelist'] # write sww intro and grid stuff. if basename_out is None: swwname = basename_in + '.sww' else: swwname = basename_out + '.sww' if verbose: print 'Output to ', swwname outfile = NetCDFFile(swwname, 'w') # For a different way of doing this, check out tsh2sww # work out sww_times and the index range this covers sww = Write_sww() sww.store_header(outfile, times, len(volumes), len(points_utm), verbose=verbose) outfile.mean_stage = mean_stage outfile.zscale = zscale sww.store_triangulation(outfile, points_utm, volumes, elevation, zone, new_origin=origin, verbose=verbose) if verbose: print 'Converting quantities' j = 0 # Read in a time slice from each mux file and write it to the sww file for ha, ua, va in map(None, mux['HA'], mux['UA'], mux['VA']): if j >= mux_times_start_i and j < mux_times_fin_i: stage = zscale*ha + mean_stage h = stage - elevation xmomentum = ua*h ymomentum = -1*va*h # -1 since in mux files south is positive. sww.store_quantities(outfile, slice_index=j - mux_times_start_i, verbose=verbose, stage=stage, xmomentum=xmomentum, ymomentum=ymomentum) j += 1 if verbose: sww.verbose_quantities(outfile) outfile.close() # # Do some conversions while writing the sww file def mux2sww_time(mux_times, mint, maxt): """ """ if mint == None: mux_times_start_i = 0 else: mux_times_start_i = searchsorted(mux_times, mint) if maxt == None: mux_times_fin_i = len(mux_times) else: maxt += 0.5 # so if you specify a time where there is # data that time will be included mux_times_fin_i = searchsorted(mux_times, maxt) return mux_times_start_i, mux_times_fin_i class Write_sww: from anuga.shallow_water.shallow_water_domain import Domain # FIXME (Ole): Hardwiring the conserved quantities like # this could be a problem. I would prefer taking them from # the instantiation of Domain. # # (DSG) There is not always a Domain instance when Write_sww is used. # Check to see if this is the same level of hardwiring as is in # shallow water doamain. sww_quantities = Domain.conserved_quantities RANGE = '_range' EXTREMA = ':extrema' def __init__(self): pass def store_header(self, outfile, times, number_of_volumes, number_of_points, description='Converted from XXX', smoothing=True, order=1, verbose=False): """ outfile - the name of the file that will be written times - A list of the time slice times OR a start time Note, if a list is given the info will be made relative. number_of_volumes - the number of triangles """ outfile.institution = 'Geoscience Australia' outfile.description = description # For sww compatibility if smoothing is True: # Smoothing to be depreciated outfile.smoothing = 'Yes' outfile.vertices_are_stored_uniquely = 'False' else: # Smoothing to be depreciated outfile.smoothing = 'No' outfile.vertices_are_stored_uniquely = 'True' outfile.order = order try: revision_number = get_revision_number() except: revision_number = None # Allow None to be stored as a string outfile.revision_number = str(revision_number) # times - A list or array of the time slice times OR a start time # times = ensure_numeric(times) # Start time in seconds since the epoch (midnight 1/1/1970) # This is being used to seperate one number from a list. # what it is actually doing is sorting lists from numeric arrays. if type(times) is list or type(times) is ArrayType: number_of_times = len(times) times = ensure_numeric(times) if number_of_times == 0: starttime = 0 else: starttime = times[0] times = times - starttime #Store relative times else: number_of_times = 0 starttime = times #times = ensure_numeric([]) outfile.starttime = starttime # dimension definitions outfile.createDimension('number_of_volumes', number_of_volumes) outfile.createDimension('number_of_vertices', 3) outfile.createDimension('numbers_in_range', 2) if smoothing is True: outfile.createDimension('number_of_points', number_of_points) # FIXME(Ole): This will cause sww files for paralle domains to # have ghost nodes stored (but not used by triangles). # To clean this up, we have to change get_vertex_values and # friends in quantity.py (but I can't be bothered right now) else: outfile.createDimension('number_of_points', 3*number_of_volumes) outfile.createDimension('number_of_timesteps', number_of_times) # variable definitions outfile.createVariable('x', precision, ('number_of_points',)) outfile.createVariable('y', precision, ('number_of_points',)) outfile.createVariable('elevation', precision, ('number_of_points',)) q = 'elevation' outfile.createVariable(q+Write_sww.RANGE, precision, ('numbers_in_range',)) # Initialise ranges with small and large sentinels. # If this was in pure Python we could have used None sensibly outfile.variables[q+Write_sww.RANGE][0] = max_float # Min outfile.variables[q+Write_sww.RANGE][1] = -max_float # Max # 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',)) for q in Write_sww.sww_quantities: outfile.createVariable(q, precision, ('number_of_timesteps', 'number_of_points')) outfile.createVariable(q+Write_sww.RANGE, precision, ('numbers_in_range',)) # Initialise ranges with small and large sentinels. # If this was in pure Python we could have used None sensibly outfile.variables[q+Write_sww.RANGE][0] = max_float # Min outfile.variables[q+Write_sww.RANGE][1] = -max_float # Max if type(times) is list or type(times) is ArrayType: outfile.variables['time'][:] = times #Store time relative if verbose: print '------------------------------------------------' print 'Statistics:' print ' t in [%f, %f], len(t) == %d'\ %(min(times.flat), max(times.flat), len(times.flat)) def store_triangulation(self, outfile, points_utm, volumes, elevation, zone=None, new_origin=None, points_georeference=None, verbose=False): """ new_origin - qa georeference that the points can be set to. (Maybe do this before calling this function.) points_utm - currently a list or array of the points in UTM. points_georeference - the georeference of the points_utm How about passing new_origin and current_origin. If you get both, do a convertion from the old to the new. If you only get new_origin, the points are absolute, convert to relative if you only get the current_origin the points are relative, store as relative. if you get no georefs create a new georef based on the minimums of points_utm. (Another option would be to default to absolute) Yes, and this is done in another part of the code. Probably geospatial. If you don't supply either geo_refs, then supply a zone. If not the default zone will be used. precon header has been called. """ number_of_points = len(points_utm) volumes = array(volumes) points_utm = array(points_utm) # given the two geo_refs and the points, do the stuff # described in the method header # if this is needed else where, pull out as a function points_georeference = ensure_geo_reference(points_georeference) new_origin = ensure_geo_reference(new_origin) if new_origin is None and points_georeference is not None: points = points_utm geo_ref = points_georeference else: if new_origin is None: new_origin = Geo_reference(zone,min(points_utm[:,0]), min(points_utm[:,1])) points = new_origin.change_points_geo_ref(points_utm, points_georeference) geo_ref = new_origin # At this stage I need a georef and points # the points are relative to the georef geo_ref.write_NetCDF(outfile) # This will put the geo ref in the middle #geo_ref=Geo_reference(refzone,(max(x)+min(x))/2.0,(max(x)+min(y))/2.) x = points[:,0] y = points[:,1] z = outfile.variables['z'][:] if verbose: print '------------------------------------------------' print 'More Statistics:' print ' Extent (/lon):' print ' x in [%f, %f], len(lat) == %d'\ %(min(x), max(x), len(x)) print ' y in [%f, %f], len(lon) == %d'\ %(min(y), max(y), len(y)) print ' z in [%f, %f], len(z) == %d'\ %(min(elevation), max(elevation), len(elevation)) print 'geo_ref: ',geo_ref print '------------------------------------------------' #z = resize(bath_grid,outfile.variables['z'][:].shape) outfile.variables['x'][:] = points[:,0] #- geo_ref.get_xllcorner() outfile.variables['y'][:] = points[:,1] #- geo_ref.get_yllcorner() outfile.variables['z'][:] = elevation outfile.variables['elevation'][:] = elevation #FIXME HACK outfile.variables['volumes'][:] = volumes.astype(Int32) #On Opteron 64 q = 'elevation' # This updates the _range values outfile.variables[q+Write_sww.RANGE][0] = min(elevation) outfile.variables[q+Write_sww.RANGE][1] = max(elevation) def store_quantities(self, outfile, precision=Float, slice_index=None, time=None, verbose=False, **quant): """ Write the quantity info. **quant is extra keyword arguments passed in. These must be the sww quantities, currently; stage, xmomentum, ymomentum. if the time array is already been built, use the slice_index to specify the index. Otherwise, use time to increase the time dimension Maybe make this general, but the viewer assumes these quantities, so maybe we don't want it general - unless the viewer is general precon triangulation and header have been called. """ if time is not None: file_time = outfile.variables['time'] slice_index = len(file_time) file_time[slice_index] = time # write the conserved quantities from Domain. # Typically stage, xmomentum, ymomentum # other quantities will be ignored, silently. for q in Write_sww.sww_quantities: if not quant.has_key(q): msg = 'SWW file can not write quantity %s' %q raise NewQuantity, msg else: q_values = quant[q] outfile.variables[q][slice_index] = q_values.astype(precision) # This updates the _range values q_range = outfile.variables[q+Write_sww.RANGE][:] q_values_min = min(q_values) if q_values_min < q_range[0]: outfile.variables[q+Write_sww.RANGE][0] = q_values_min q_values_max = max(q_values) if q_values_max > q_range[1]: outfile.variables[q+Write_sww.RANGE][1] = q_values_max def verbose_quantities(self, outfile): print '------------------------------------------------' print 'More Statistics:' for q in Write_sww.sww_quantities: print ' %s in [%f, %f]' %(q, outfile.variables[q+Write_sww.RANGE][0], outfile.variables[q+Write_sww.RANGE][1]) print '------------------------------------------------' def obsolete_write_sww_time_slices(outfile, has, uas, vas, elevation, mean_stage=0, zscale=1, verbose=False): #Time stepping stage = outfile.variables['stage'] xmomentum = outfile.variables['xmomentum'] ymomentum = outfile.variables['ymomentum'] n = len(has) j=0 for ha, ua, va in map(None, has, uas, vas): if verbose and j%((n+10)/10)==0: print ' Doing %d of %d' %(j, n) w = zscale*ha + mean_stage stage[j] = w h = w - elevation xmomentum[j] = ua*h ymomentum[j] = -1*va*h # -1 since in mux files south is positive. j += 1 def urs2txt(basename_in, location_index=None): """ Not finished or tested """ files_in = [basename_in + WAVEHEIGHT_MUX_LABEL, basename_in + EAST_VELOCITY_LABEL, basename_in + NORTH_VELOCITY_LABEL] quantities = ['HA','UA','VA'] d = "," # instanciate urs_points of the three mux files. mux = {} for quantity, file in map(None, quantities, files_in): mux[quantity] = Urs_points(file) # Could check that the depth is the same. (hashing) # handle to a mux file to do depth stuff a_mux = mux[quantities[0]] # Convert to utm latitudes = a_mux.lonlatdep[:,1] longitudes = a_mux.lonlatdep[:,0] points_utm, zone = convert_from_latlon_to_utm( \ latitudes=latitudes, longitudes=longitudes) #print "points_utm", points_utm #print "zone", zone depths = a_mux.lonlatdep[:,2] # fid = open(basename_in + '.txt', 'w') fid.write("zone: " + str(zone) + "\n") if location_index is not None: #Title li = location_index fid.write('location_index'+d+'lat'+d+ 'long' +d+ 'Easting' +d+ \ 'Northing' + "\n") fid.write(str(li) +d+ str(latitudes[li])+d+ \ str(longitudes[li]) +d+ str(points_utm[li][0]) +d+ \ str(points_utm[li][01]) + "\n") # the non-time dependent stuff #Title fid.write('location_index'+d+'lat'+d+ 'long' +d+ 'Easting' +d+ \ 'Northing' +d+ 'depth m' + "\n") i = 0 for depth, point_utm, lat, long in map(None, depths, points_utm, latitudes, longitudes): fid.write(str(i) +d+ str(lat)+d+ str(long) +d+ str(point_utm[0]) +d+ \ str(point_utm[01]) +d+ str(depth) + "\n") i +=1 #Time dependent if location_index is not None: time_step = a_mux.time_step i = 0 #Title fid.write('time' +d+ 'HA depth m'+d+ \ 'UA momentum East x m/sec' +d+ 'VA momentum North y m/sec' \ + "\n") for HA, UA, VA in map(None, mux['HA'], mux['UA'], mux['VA']): fid.write(str(i*time_step) +d+ str(HA[location_index])+d+ \ str(UA[location_index]) +d+ str(VA[location_index]) \ + "\n") i +=1 class Urs_points: """ Read the info in URS mux files. for the quantities heres a correlation between the file names and what they mean; z-mux is height above sea level, m e-mux is velocity is Eastern direction, m/s n-mux is velocity is Northern direction, m/s """ def __init__(self,urs_file): self.iterated = False columns = 3 # long, lat , depth mux_file = open(urs_file, 'rb') # Number of points/stations (self.points_num,)= unpack('i',mux_file.read(4)) # nt, int - Number of time steps (self.time_step_count,)= unpack('i',mux_file.read(4)) #print "self.time_step_count", self.time_step_count #dt, float - time step, seconds (self.time_step,) = unpack('f', mux_file.read(4)) #print "self.time_step", self.time_step msg = "Bad data in the urs file." if self.points_num < 0: mux_file.close() raise ANUGAError, msg if self.time_step_count < 0: mux_file.close() raise ANUGAError, msg if self.time_step < 0: mux_file.close() raise ANUGAError, msg # the depth is in meters, and it is the distance from the ocean # to the sea bottom. lonlatdep = p_array.array('f') lonlatdep.read(mux_file, columns * self.points_num) lonlatdep = array(lonlatdep, typecode=Float) lonlatdep = reshape(lonlatdep, (self.points_num, columns)) #print 'lonlatdep',lonlatdep self.lonlatdep = lonlatdep self.mux_file = mux_file # check this array def __iter__(self): """ iterate over quantity data which is with respect to time. Note: You can only interate once over an object returns quantity infomation for each time slice """ msg = "You can only interate once over a urs file." assert not self.iterated, msg self.iter_time_step = 0 self.iterated = True return self def next(self): if self.time_step_count == self.iter_time_step: self.close() raise StopIteration #Read in a time slice from mux file hz_p_array = p_array.array('f') hz_p_array.read(self.mux_file, self.points_num) hz_p = array(hz_p_array, typecode=Float) self.iter_time_step += 1 return hz_p def close(self): self.mux_file.close() #### END URS UNGRIDDED 2 SWW ### def start_screen_catcher(dir_name, myid='', numprocs='', extra_info='', print_to_screen=False, verbose=False): """ Used to store screen output and errors to file, if run on multiple processes eachprocessor will have its own output and error file. extra_info - is used as a string that can identify outputs with another string eg. '_other' """ import sys dir_name = dir_name if access(dir_name,W_OK) == 0: if verbose: print 'Make directory %s' %dir_name if verbose: print "myid", myid mkdir (dir_name,0777) if myid <>'': myid = '_'+str(myid) if numprocs <>'': numprocs = '_'+str(numprocs) if extra_info <>'': extra_info = '_'+str(extra_info) screen_output_name = dir_name + "screen_output%s%s%s.txt" %(myid, numprocs, extra_info) screen_error_name = dir_name + "screen_error%s%s%s.txt" %(myid, numprocs, extra_info) print screen_output_name #used to catch screen output to file sys.stdout = Screen_Catcher(screen_output_name) sys.stderr = Screen_Catcher(screen_error_name) class Screen_Catcher: """this simply catches the screen output and stores it to file defined by start_screen_catcher (above) """ def __init__(self, filename): self.filename = filename if exists(self.filename)is True: print'Old existing file "%s" has been deleted' %(self.filename) remove(self.filename) def write(self, stuff): fid = open(self.filename, 'a') fid.write(stuff) # if print_to_screen: print stuff def copy_code_files(dir_name, filename1, filename2): """Copies "filename1" and "filename2" to "dir_name". Very useful for information management filename1 and filename2 are both absolute pathnames """ if access(dir_name,F_OK) == 0: print 'Make directory %s' %dir_name mkdir (dir_name,0777) shutil.copy(filename1, dir_name + sep + basename(filename1)) shutil.copy(filename2, dir_name + sep + basename(filename2)) # copy (__file__, project.output_run_time_dir + basename(__file__)) print 'Files %s and %s copied' %(filename1, filename2) def get_data_from_file(filename,separator_value = ','): """ Read in data information from file and Returns: header_fields, a string? of the first line separated by the 'separator_value' data, a array (N data columns X M lines) in the file excluding the header NOTE: wont deal with columns with different lenghts and there must be no blank lines at the end. """ fid = open(filename) lines = fid.readlines() fid.close() header_line = lines[0] header_fields = header_line.split(separator_value) #array to store data, number in there is to allow float... #i'm sure there is a better way! data=array([],typecode=Float) data=resize(data,((len(lines)-1),len(header_fields))) # print 'number of fields',range(len(header_fields)) # print 'number of lines',len(lines), shape(data) # print'data',data[1,1],header_line array_number = 0 line_number = 1 while line_number < (len(lines)): for i in range(len(header_fields)): #this get line below the header, explaining the +1 #and also the line_number can be used as the array index fields = lines[line_number].split(separator_value) #assign to array data[array_number,i] = float(fields[i]) line_number = line_number +1 array_number = array_number +1 return header_fields, data def store_parameters(verbose=False,**kwargs): """ Store "kwargs" into a temp csv file, if "completed" is a kwargs csv file is kwargs[file_name] else it is kwargs[output_dir] + details_temp.csv Must have a file_name keyword arg, this is what is writing to. might be a better way to do this using CSV module Writer and writeDict writes file to "output_dir" unless "completed" is in kwargs, then it writes to "file_name" kwargs """ import types # import os # Check that kwargs is a dictionary if type(kwargs) != types.DictType: raise TypeError #is completed is kwargs? try: kwargs['completed'] completed=True except: completed=False #get file name and removes from dict and assert that a file_name exists if completed: try: file = str(kwargs['file_name']) except: raise 'kwargs must have file_name' else: #write temp file in output directory try: file = str(kwargs['output_dir'])+'detail_temp.csv' except: raise 'kwargs must have output_dir' #extracts the header info and the new line info line='' header='' count=0 keys = kwargs.keys() keys.sort() #used the sorted keys to create the header and line data for k in keys: # print "%s = %s" %(k, kwargs[k]) header = header+str(k) line = line+str(kwargs[k]) count+=1 if count 0 Usage: max_runup = get_maximum_inundation_elevation(filename, polygon=None, time_interval=None, verbose=False) filename is a NetCDF sww file containing ANUGA model output. Optional arguments polygon and time_interval restricts the maximum runup calculation to a points that lie within the specified polygon and time interval. If no inundation is found within polygon and time_interval the return value is None signifying "No Runup" or "Everything is dry". See general function get_maximum_inundation_data for details. """ runup, _ = get_maximum_inundation_data(filename, polygon=polygon, time_interval=time_interval, verbose=verbose) return runup def get_maximum_inundation_location(filename, polygon=None, time_interval=None, verbose=False): """Return location of highest elevation where h > 0 Usage: max_runup_location = get_maximum_inundation_location(filename, polygon=None, time_interval=None, verbose=False) filename is a NetCDF sww file containing ANUGA model output. Optional arguments polygon and time_interval restricts the maximum runup calculation to a points that lie within the specified polygon and time interval. If no inundation is found within polygon and time_interval the return value is None signifying "No Runup" or "Everything is dry". See general function get_maximum_inundation_data for details. """ _, max_loc = get_maximum_inundation_data(filename, polygon=polygon, time_interval=time_interval, verbose=verbose) return max_loc def get_maximum_inundation_data(filename, polygon=None, time_interval=None, use_centroid_values=False, verbose=False): """Compute maximum run up height from sww file. Usage: runup, location = get_maximum_inundation_data(filename, polygon=None, time_interval=None, verbose=False) Algorithm is as in get_maximum_inundation_elevation from shallow_water_domain except that this function works with the sww file and computes the maximal runup height over multiple timesteps. Optional arguments polygon and time_interval restricts the maximum runup calculation to a points that lie within the specified polygon and time interval. Polygon is assumed to be in (absolute) UTM coordinates in the same zone as domain. If no inundation is found within polygon and time_interval the return value is None signifying "No Runup" or "Everything is dry". """ # We are using nodal values here as that is what is stored in sww files. # Water depth below which it is considered to be 0 in the model # FIXME (Ole): Allow this to be specified as a keyword argument as well from anuga.utilities.polygon import inside_polygon from anuga.config import minimum_allowed_height from Scientific.IO.NetCDF import NetCDFFile dir, base = os.path.split(filename) iterate_over = get_all_swwfiles(dir,base) # Read sww file if verbose: print 'Reading from %s' %filename # FIXME: Use general swwstats (when done) maximal_runup = None maximal_runup_location = None for file, swwfile in enumerate (iterate_over): # Read sww file filename = join(dir,swwfile+'.sww') if verbose: print 'Reading from %s' %filename # FIXME: Use general swwstats (when done) fid = NetCDFFile(filename) # Get geo_reference # sww files don't have to have a geo_ref try: geo_reference = Geo_reference(NetCDFObject=fid) except AttributeError, e: geo_reference = Geo_reference() # Default georef object xllcorner = geo_reference.get_xllcorner() yllcorner = geo_reference.get_yllcorner() zone = geo_reference.get_zone() # Get extent volumes = fid.variables['volumes'][:] x = fid.variables['x'][:] + xllcorner y = fid.variables['y'][:] + yllcorner # Get the relevant quantities elevation = fid.variables['elevation'][:] stage = fid.variables['stage'][:] # Here's where one could convert nodal information to centroid # information # but is probably something we need to write in C. # Here's a Python thought which is NOT finished!!! if use_centroid_values is True: x = get_centroid_values(x, volumes) y = get_centroid_values(y, volumes) elevation = get_centroid_values(elevation, volumes) # Spatial restriction if polygon is not None: msg = 'polygon must be a sequence of points.' assert len(polygon[0]) == 2, msg # FIXME (Ole): Make a generic polygon input check in polygon.py # and call it here points = concatenate((x[:,NewAxis], y[:,NewAxis]), axis=1) point_indices = inside_polygon(points, polygon) # Restrict quantities to polygon elevation = take(elevation, point_indices) stage = take(stage, point_indices, axis=1) # Get info for location of maximal runup points_in_polygon = take(points, point_indices) x = points_in_polygon[:,0] y = points_in_polygon[:,1] else: # Take all points point_indices = arange(len(x)) # Temporal restriction time = fid.variables['time'][:] all_timeindices = arange(len(time)) if time_interval is not None: msg = 'time_interval must be a sequence of length 2.' assert len(time_interval) == 2, msg msg = 'time_interval %s must not be decreasing.' %(time_interval) assert time_interval[1] >= time_interval[0], msg msg = 'Specified time interval [%.8f:%.8f]' %tuple(time_interval) msg += ' must does not match model time interval: [%.8f, %.8f]\n'\ %(time[0], time[-1]) if time_interval[1] < time[0]: raise ValueError(msg) if time_interval[0] > time[-1]: raise ValueError(msg) # Take time indices corresponding to interval (& is bitwise AND) timesteps = compress((time_interval[0] <= time) & (time <= time_interval[1]), all_timeindices) msg = 'time_interval %s did not include any model timesteps.' %(time_interval) assert not alltrue(timesteps == 0), msg else: # Take them all timesteps = all_timeindices fid.close() # Compute maximal runup for each timestep #maximal_runup = None #maximal_runup_location = None #maximal_runups = [None] #maximal_runup_locations = [None] for i in timesteps: if use_centroid_values is True: stage_i = get_centroid_values(stage[i,:], volumes) else: stage_i = stage[i,:] depth = stage_i - elevation # Get wet nodes i.e. nodes with depth>0 within given region and timesteps wet_nodes = compress(depth > minimum_allowed_height, arange(len(depth))) if alltrue(wet_nodes == 0): runup = None else: # Find maximum elevation among wet nodes wet_elevation = take(elevation, wet_nodes) runup_index = argmax(wet_elevation) runup = max(wet_elevation) assert wet_elevation[runup_index] == runup # Must always be True if runup > maximal_runup: maximal_runup = runup # This works even if maximal_runups is None #print "NEW RUNUP",runup # Record location wet_x = take(x, wet_nodes) wet_y = take(y, wet_nodes) maximal_runup_location = [wet_x[runup_index], wet_y[runup_index]] #print 'maximal_runup, maximal_runup_location',maximal_runup, maximal_runup_location return maximal_runup, maximal_runup_location def get_all_swwfiles(look_in_dir='',base_name='',verbose=False): ''' Finds all the sww files in a "look_in_dir" which contains a "base_name". will accept base_name with or without the extension ".sww" Returns: a list of strings Usage: iterate_over = get_all_swwfiles(dir, name) then for swwfile in iterate_over: do stuff Check "export_grids" and "get_maximum_inundation_data" for examples ''' #plus tests the extension name, extension = os.path.splitext(base_name) if extension <>'' and extension <> '.sww': msg = msg='file %s %s must be an NetCDF sww file!'%(base_name,extension) raise IOError, msg if look_in_dir == "": look_in_dir = "." # Unix compatibility dir_ls = os.listdir(look_in_dir) #print 'dir_ls',dir_ls, base iterate_over = [x[:-4] for x in dir_ls if name in x and x[-4:] == '.sww'] if len(iterate_over) == 0: msg = 'No files of the base name %s'\ %(name) raise IOError, msg if verbose: print 'iterate over %s' %(iterate_over) return iterate_over def get_all_files_with_extension(look_in_dir='',base_name='',extension='.sww',verbose=False): ''' Finds all the sww files in a "look_in_dir" which contains a "base_name". Returns: a list of strings Usage: iterate_over = get_all_swwfiles(dir, name) then for swwfile in iterate_over: do stuff Check "export_grids" and "get_maximum_inundation_data" for examples ''' #plus tests the extension name, ext = os.path.splitext(base_name) # print 'look_in_dir',look_in_dir if ext <>'' and ext <> extension: msg = msg='base_name %s must be an file with %s extension!'%(base_name,extension) raise IOError, msg if look_in_dir == "": look_in_dir = "." # Unix compatibility # print 'look_in_dir',look_in_dir, getcwd() dir_ls = os.listdir(look_in_dir) #print 'dir_ls',dir_ls, base_name iterate_over = [x[:-4] for x in dir_ls if name in x and x[-4:] == extension] if len(iterate_over) == 0: msg = 'No files of the base name %s in %s'\ %(name, look_in_dir) raise IOError, msg if verbose: print 'iterate over %s' %(iterate_over) return iterate_over def get_all_directories_with_name(look_in_dir='',base_name='',verbose=False): ''' Finds all the sww files in a "look_in_dir" which contains a "base_name". Returns: a list of strings Usage: iterate_over = get_all_swwfiles(dir, name) then for swwfile in iterate_over: do stuff Check "export_grids" and "get_maximum_inundation_data" for examples ''' #plus tests the extension if look_in_dir == "": look_in_dir = "." # Unix compatibility # print 'look_in_dir',look_in_dir dir_ls = os.listdir(look_in_dir) # print 'dir_ls',dir_ls iterate_over = [x for x in dir_ls if base_name in x] if len(iterate_over) == 0: msg = 'No files of the base name %s'\ %(name) raise IOError, msg if verbose: print 'iterate over %s' %(iterate_over) return iterate_over #------------------------------------------------------------- if __name__ == "__main__": #setting umask from config to force permissions for all files and directories # created to the same. (it was noticed the "mpirun" doesn't honour the umask # set in your .bashrc etc file) from config import umask import os os.umask(umask)