""" --------------------------------------------------------------get_tuflow_data.py ANUGA UTILITIES TO READ TUFLOW MID/MIF INPUT FILES These functions assist Anuga to read in data from a previous tuflow model get_tuflow_bc Reads in the inflow and dsbc timeseries get_wbnm_coords Gets the outlet coords of each inflow subarea get_tuflow_bc_filetype Gets the tuflow (ts1 or csv) bc filetype convert_asc_ts_to_tms Converts an ascii timeseries to NetCDF format get_tuflow_tmf Reads in the surface classes and properties get_tuflow_mat Reads in the surface class polygons get_tuflow_2d_po Reads the tuflow po flowlines This script is imported by model_run.py Code Version : 1.00 June 2008 Initial release Author : E Rigby (0437 250 500) ted.rigby@rienco.com.au ------------------------------------------------------------------------------- """ import string #START FUNCTION ================================================= get_tuflow_bc def get_tuflow_bc(location,loc_filename,bc_filetype,verbose): 'returns an ascii hydrograph for that location' if verbose : print " ++++ get_tuflow_bc - Opening file %s to obtain temporal data for %s type bc at %s " %(loc_filename,bc_filetype,location) if bc_filetype=="ts1": # Is a ts1 florate file - open and read in the event and location specific flow hydrograph bctype = "inflow" if verbose : print " ++++ get_tuflow_bc - bc type is %s" %(bctype) ts1_fid=open(loc_filename, 'r') # open the associated hydrograph file for that loc ts1_lines=ts1_fid.readlines() # read all file into string list ts1_fid.close() # close the loc specific hydrograph file subsline=ts1_lines[4] # select the line with subarea headings if verbose : print " ++++ get_tuflow_bc - Subarea header line is ", subsline [:50] index=subsline.index(location)+len(location) -12 # locate start pos of selected subarea header col tsecs=[] qmecs=[] for ts1_line in ts1_lines[5:]: t=float(ts1_line[:12])*60.0 # time convert time to secs - was minutes if (ts1_line[index:index+12]) == " ": q=0 # end of this hydrograph not all full duration else: q=float(ts1_line[index:index+12])# flowrate in m3/sec tsecs.append(t) qmecs.append(q) if verbose : print " ++++ get_tuflow_bc - ts1 flow hydrograph at %s (first 5 values)" %(location) for i in range (5): print " ++++ ", tsecs[i],qmecs[i] return tsecs,qmecs else: # Is a dsbc csv file - open and read in the event and location specific elevation hydrograph bctype = "dsbc" if verbose : print " ++++ bc type is %s" %(bctype) csv_fid=open(loc_filename, 'r') # open the associated elev hydrograph file for that loc csv_lines=csv_fid.readlines() # read all file into string list csv_fid.close() # close the loc specific hydrograph file tsecs=[] elev =[] for csv_line in csv_lines[1:]: # read the dsbc data sep_col=csv_line.find(',') # locate seperator - not fixed field widths t=float(csv_line[:sep_col])* 3600.0 # time convert to seconds for Anuga - was hours e=float(csv_line[sep_col+1:])# elevation at bdry in (m) tsecs.append(t) elev.append(e) if verbose : for i in range (5): print " ++++ csv elev hydrograph at ", location,tsecs[i],elev[i] return tsecs,elev #END FUNCTION =================================================== get_tuflow_bc #START FUNCTION =============================================== get_wbnm_coords def get_wbnm_coords(wbnm_filename,subarea,verbose): 'reurns the x and y outlet coordinates of the subarea' wbn_fid=open(wbnm_filename, 'r') # open the wbnm run file for wbn_line in wbn_fid: # read file a line at a time if wbn_line[:12].find(subarea.strip())>=0: # found the subarea line xcoord = float(wbn_line[36:47]) # Subareas outlet xcoord ycoord = float(wbn_line[48:59]) # Subareas outlet xcoord wbn_fid.close() # close the wbn run file if verbose : print ' ++++ get_wbnm_coords - Subareal outlet for %s found at %7.2f %7.2f ' %(subarea,xcoord,ycoord) return xcoord, ycoord # return the outlet coords if wbn_line.find("#####END_TOPOLOGY")>0: # did not find the subarea in topology print " >>>> ERROR - get_wbnm_coords - Could not find %s in %s" %(subarea,wbnm_filename), sys.exit(0) # terminate the program #END FUNCTION =================================================== get_wbnm_coords #START FUNCTION =========================================== get_tuflow_bc_filetype def get_tuflow_bc_filetype(loc_filename,verbose): 'returns type ("ts1") ("csv") of bc_file' if loc_filename.find(".ts1")>=0: # Is file containing hydrographs for each subarea bc_filetype="ts1" return bc_filetype elif loc_filename.find(".csv")>=0: # Is time(hrs),elev(m) dsbc.csv file bc_filetype="csv" return bc_filetype else: print " >>>> ERROR get_tuflow_bc_filetype - Unrecognised bc filetype %s - aborting " %(bc_filetype) sys.exit(0) # terminate the program #END FUNCTION =========================================== get_tuflow_bc_filetype #START FUNCTION ========================================== convert_asc_ts_to_tms def convert_asc_ts_to_tms(tms_dir,ts_name, time_units, time, value,verbose): 'saves the [time][value] timeseries into a netcdf tms object ' import os from os import sep from Numeric import array, Float from Scientific.IO.NetCDF import NetCDFFile if verbose : print " ++++ asc_ts_to_tms - Converting %s time series in %s to tms " %(ts_name,time_units) for i in range (5): # len(time) would be full list print " ++++ ", time[i],value[i] tms_filename = tms_dir + sep + ts_name.strip() + '.tms' # buildname of tms file to be created if verbose : print " ++++ asc_ts_to_tms - Created tms file named %s to store time series in" %(tms_filename) # compute time conversion factor to bring to secs if time_units=='days': factor=60.0*60.0*24.0 elif time_units=='hours': factor=60.0*60.0 elif time_units=='minutes': factor=60.0 else: factor=1.0 for i in range (len(time)): time[i]=time[i]*factor # convert to and create NetCDF tms file of time series T=array(time,Float) # secs V=array(value,Float) # whatever tms_fid=NetCDFFile(tms_filename,'w') tms_fid.institution='Rienco Consulting' tms_fid.description='Utility to create NetCDF (tms) file of ascii time series' tms_fid.starttime=0.0 tms_fid.createDimension('number_of_timesteps', len(T)) tms_fid.createVariable('time',Float,('number_of_timesteps',) ) tms_fid.variables['time'][:] = T tms_fid.createVariable('value',Float,('number_of_timesteps',) ) tms_fid.variables['value'][:]= V tms_fid.close() if verbose : print " ++++ asc_ts_to_tms - NetCDF time series created and stored in %s - time units are %s " %(tms_filename,time_units) return tms_filename #END FUNCTION ========================================== convert_asc_ts_to_tms #START FUNCTION =============================================== get_tuflow_tmf def get_tuflow_tmf(tmf_filename,verbose): 'read in the tuflow roughness values and return as python list of values' tmf_fid=open(tmf_filename, 'r') # open the tuflow material (roughness) file tmf_lines=tmf_fid.readlines() # read all file into string list tmf_fid.close() # close the tmf file mat_roughness_data_list=[] if verbose : print " ++++ get_tuflow_tmf - Surface roughness data read in from tmf and stored in mat_roughness_data_list " for tmf_line in tmf_lines: #if tmf_line[0]=='!' or tmf_line[0]==' ' or tmf_line[0]=='\n': if tmf_line[0]=='!' or tmf_line[0] in string.whitespace: pass else: # decode and store the list values tmf_line=tmf_line.replace('!',',') # get rid of the ! so can use split tmf_line_fields = tmf_line.split(',') # break the line into fields tmf_line_fields[0] = int(tmf_line_fields[0]) # The tuflow surface RID tmf_line_fields[1] = float(tmf_line_fields[1]) # The tuflow fixed n value tmf_line_fields[2] = int(tmf_line_fields[2]) # Not used tmf_line_fields[3] = int(tmf_line_fields[3]) # Not used tmf_line_fields[4] = float(tmf_line_fields[4]) # The tuflow d1 value tmf_line_fields[5] = float(tmf_line_fields[5]) # The tuflow n1 value tmf_line_fields[6] = float(tmf_line_fields[6]) # The tuflow d2 value tmf_line_fields[7] = float(tmf_line_fields[7]) # The tuflow n2 value tmf_line_fields[8] = str(tmf_line_fields[8]) # The tuflow surface description if verbose : print " ++++ get_tuflow_tmf appending " ,tmf_line_fields mat_roughness_data_list.append(tmf_line_fields) # there are nine data fields including the description return mat_roughness_data_list # a list of surface data fields as in the tuflow tmf file #END FUNCTION ================================================ get_tuflow_tmf #START FUNCTION =============================================== get_tuflow_mat def get_tuflow_mat(mat_filename,verbose): 'reads in the roughness polys from the tuflow 2d_mat.mif file and IDs from the mid' # Read in the polygons from the 2d_mat.mif file matmif_fid=open(mat_filename, 'r') # open the tuflow material (roughness polys) file matmif_lines=matmif_fid.readlines() # read all file into string list matmif_fid.close() # close the 2d_mat.mif file poly_coords_list=[] # float list of roughness polygon coord lists point_coords=[] # float list of single coord x y pair for i in range (len(matmif_lines)): if matmif_lines[i][:12].find('Region')>=0: nlines=int(matmif_lines[i+1] ) # read next line get no of coords (nlines) in poly poly=[] # float list of coords x,y for j in range (nlines): # read next nlines of xy coords into poly[] x,y = matmif_lines[i+2+j].split() # split into two strings x = float(x) # convert to numeric coords y = float(y) point_coords = [x,y] poly.append(point_coords) # add to poly as float xy pairs poly_coords_list.append(poly) # append this poly coord list to list of polys i=i+1+j # reset counter to previous to last line read to avoid problem when at end of file else: # unwanted - skip line pass # Read in the polygon material (roughness) IDs from the 2d_mat.mid file matmid_fid=open(mat_filename[:-4] + '.mid', 'r') # open the tuflow material (roughness polys) file matmid_lines=matmid_fid.readlines() # read all file into string list matmid_fid.close() # close the 2d_mat.mid file poly_RID_list=[] # integer list identifying roughness polygons ID for matmid_line in matmid_lines: RID=int(matmid_line[:(matmid_line.index(','))]) # get the RID of each poly poly_RID_list.append(RID) if len(poly_RID_list)== len(poly_coords_list): # all as it should be if verbose : print " ++++ get_tuflow_mat - Read in %i surface roughness polygons from %s " %(len(poly_RID_list), mat_filename) return poly_RID_list,poly_coords_list # [intIDs],[list of poly floatcoords list] else: print " >>>> ERROR get_tuflow_mat - Poly regions list and poly ID list not of same length in get_tuflow_mat??" print " >>>> ERROR get_tuflow_mat - Poly ID list %d and poly list %d" %(len(mat_ID_list),len(mat_poly_list)) sys.exit(0) # terminate the program #END FUNCTION ================================================ get_tuflow_mat