source: anuga_work/development/ted_rigby/get_tuflow_data.py @ 5977

Last change on this file since 5977 was 5973, checked in by ole, 15 years ago

Ted Rigby's variable friction code

File size: 13.4 KB
Line 
1"""
2--------------------------------------------------------------get_tuflow_data.py
3ANUGA UTILITIES TO READ TUFLOW MID/MIF INPUT FILES
4These functions assist Anuga to read in data from a previous tuflow model
5
6            get_tuflow_bc            Reads in the inflow and dsbc timeseries
7            get_wbnm_coords          Gets the outlet coords of each inflow subarea
8            get_tuflow_bc_filetype   Gets the tuflow (ts1 or csv) bc filetype
9            convert_asc_ts_to_tms    Converts an ascii timeseries to NetCDF format
10            get_tuflow_tmf           Reads in the  surface classes and properties
11            get_tuflow_mat           Reads in the surface class polygons
12            get_tuflow_2d_po         Reads the tuflow po flowlines
13           
14           
15
16This script is imported by model_run.py
17
18Code Version : 1.00 June 2008 Initial  release
19Author       : E Rigby (0437 250 500) ted.rigby@rienco.com.au
20
21-------------------------------------------------------------------------------
22"""
23
24import string
25
26#START FUNCTION ================================================= get_tuflow_bc
27
28def get_tuflow_bc(location,loc_filename,bc_filetype,verbose):
29    'returns an ascii hydrograph for that location'
30               
31    if verbose :
32        print "         ++++ get_tuflow_bc - Opening file %s to obtain temporal data for %s type bc at %s " %(loc_filename,bc_filetype,location) 
33
34    if bc_filetype=="ts1":
35   
36        #   Is a ts1 florate file - open and read in the event and location specific flow hydrograph
37        bctype = "inflow"
38       
39        if verbose :
40            print "         ++++ get_tuflow_bc - bc type is  %s" %(bctype)
41        ts1_fid=open(loc_filename, 'r') # open the associated hydrograph file for that loc
42        ts1_lines=ts1_fid.readlines()   # read all file into string list
43        ts1_fid.close()                 # close the loc specific hydrograph file
44        subsline=ts1_lines[4]           # select the line with subarea headings
45
46        if verbose :
47            print "         ++++ get_tuflow_bc - Subarea header line is ", subsline [:50]
48        index=subsline.index(location)+len(location) -12  # locate start pos of selected subarea header col
49        tsecs=[]
50        qmecs=[]
51        for ts1_line in ts1_lines[5:]:
52            t=float(ts1_line[:12])*60.0          # time convert time to secs - was minutes
53            if (ts1_line[index:index+12]) == "            ":
54                q=0                              # end of this hydrograph not all full duration
55            else:
56                q=float(ts1_line[index:index+12])# flowrate in m3/sec
57            tsecs.append(t)
58            qmecs.append(q)
59        if verbose :
60            print "         ++++ get_tuflow_bc - ts1 flow hydrograph at %s (first 5 values)" %(location)
61            for i in range (5):
62                print "         ++++ ", tsecs[i],qmecs[i]
63               
64        return tsecs,qmecs
65
66    else:
67
68        #    Is a dsbc csv file - open and read in the event and location specific elevation hydrograph
69        bctype = "dsbc"
70
71        if verbose :
72            print "         ++++ bc type is  %s" %(bctype)
73        csv_fid=open(loc_filename, 'r')  # open the associated elev hydrograph file for that loc
74        csv_lines=csv_fid.readlines()    # read all file into string list
75        csv_fid.close()                  # close the loc specific hydrograph file
76        tsecs=[]
77        elev =[]
78        for csv_line in csv_lines[1:]:   # read the dsbc data
79            sep_col=csv_line.find(',')   # locate seperator - not fixed field widths
80            t=float(csv_line[:sep_col])* 3600.0 # time convert to seconds for Anuga - was hours
81            e=float(csv_line[sep_col+1:])# elevation at bdry in (m)
82            tsecs.append(t)
83            elev.append(e)
84        if verbose :
85            for i in range (5):
86                print "         ++++ csv elev hydrograph at ", location,tsecs[i],elev[i]
87
88        return tsecs,elev
89
90#END FUNCTION =================================================== get_tuflow_bc       
91
92
93
94#START FUNCTION =============================================== get_wbnm_coords
95
96def get_wbnm_coords(wbnm_filename,subarea,verbose):
97    'reurns the x and y outlet coordinates of the subarea'
98   
99    wbn_fid=open(wbnm_filename, 'r')           # open the wbnm run file
100    for wbn_line in  wbn_fid:                  # read file a line at a time
101       
102        if wbn_line[:12].find(subarea.strip())>=0:     # found the subarea line
103            xcoord = float(wbn_line[36:47])    # Subareas outlet xcoord
104            ycoord = float(wbn_line[48:59])    # Subareas outlet xcoord
105            wbn_fid.close()                    # close the wbn run file
106            if verbose :
107                print '         ++++ get_wbnm_coords - Subareal outlet for %s found at %7.2f %7.2f ' %(subarea,xcoord,ycoord)   
108            return xcoord, ycoord              # return the outlet coords   
109                       
110        if wbn_line.find("#####END_TOPOLOGY")>0:   # did not find the subarea in topology
111            print "        >>>> ERROR - get_wbnm_coords - Could not find %s in %s" %(subarea,wbnm_filename),
112            sys.exit(0)                        # terminate the program
113   
114
115
116#END FUNCTION =================================================== get_wbnm_coords
117
118
119#START FUNCTION =========================================== get_tuflow_bc_filetype
120
121def get_tuflow_bc_filetype(loc_filename,verbose):
122    'returns type ("ts1") ("csv") of bc_file'
123
124    if loc_filename.find(".ts1")>=0:    # Is file containing hydrographs for each subarea
125        bc_filetype="ts1"
126        return bc_filetype
127       
128    elif loc_filename.find(".csv")>=0:  # Is time(hrs),elev(m) dsbc.csv file
129        bc_filetype="csv"
130        return bc_filetype
131       
132    else:
133        print "         >>>> ERROR get_tuflow_bc_filetype - Unrecognised bc filetype %s - aborting " %(bc_filetype)
134        sys.exit(0)                    # terminate the program
135
136
137
138#END FUNCTION =========================================== get_tuflow_bc_filetype
139
140#START FUNCTION ========================================== convert_asc_ts_to_tms
141def convert_asc_ts_to_tms(tms_dir,ts_name, time_units, time, value,verbose):
142    'saves the [time][value] timeseries into a netcdf tms object '
143
144    import os
145    from os import sep
146    from Numeric import array, Float
147    from Scientific.IO.NetCDF import NetCDFFile
148     
149    if verbose :
150        print "         ++++ asc_ts_to_tms - Converting %s time series in %s to tms " %(ts_name,time_units)
151        for i in range (5):     # len(time) would be full list
152            print "         ++++ ", time[i],value[i]
153   
154    tms_filename = tms_dir + sep + ts_name.strip() + '.tms' # buildname of tms file to be created
155    if verbose :
156        print "         ++++ asc_ts_to_tms - Created tms file named %s to store time series in" %(tms_filename)
157 
158   
159# compute time conversion factor to bring to secs
160
161    if time_units=='days':
162        factor=60.0*60.0*24.0
163    elif time_units=='hours':
164        factor=60.0*60.0
165    elif time_units=='minutes':
166        factor=60.0
167    else:
168        factor=1.0
169    for i in range (len(time)):
170        time[i]=time[i]*factor
171       
172    # convert to and create NetCDF tms file of time series
173   
174
175    T=array(time,Float)          # secs
176    V=array(value,Float)         # whatever
177
178    tms_fid=NetCDFFile(tms_filename,'w')
179    tms_fid.institution='Rienco Consulting'
180    tms_fid.description='Utility to create NetCDF (tms) file of ascii time series'
181    tms_fid.starttime=0.0
182    tms_fid.createDimension('number_of_timesteps', len(T))
183    tms_fid.createVariable('time',Float,('number_of_timesteps',) )
184    tms_fid.variables['time'][:] = T
185    tms_fid.createVariable('value',Float,('number_of_timesteps',) )
186    tms_fid.variables['value'][:]= V
187    tms_fid.close()
188   
189    if verbose :
190            print "         ++++ asc_ts_to_tms - NetCDF time series created and stored in %s - time units are %s " %(tms_filename,time_units)
191       
192    return tms_filename
193   
194
195#END FUNCTION ========================================== convert_asc_ts_to_tms
196
197#START FUNCTION =============================================== get_tuflow_tmf
198def get_tuflow_tmf(tmf_filename,verbose):
199    'read in the tuflow roughness values and return as python list of values'
200
201   
202    tmf_fid=open(tmf_filename, 'r')           # open the tuflow material (roughness) file
203    tmf_lines=tmf_fid.readlines()             # read all file into string list
204    tmf_fid.close()                           # close the tmf file
205
206    mat_roughness_data_list=[]
207
208    if verbose :
209        print "         ++++ get_tuflow_tmf - Surface roughness data read in from tmf and stored in mat_roughness_data_list " 
210
211    for tmf_line in tmf_lines:
212        #if tmf_line[0]=='!' or tmf_line[0]==' ' or tmf_line[0]=='\n':
213        if tmf_line[0]=='!' or tmf_line[0] in string.whitespace:
214            pass
215        else:                                                 # decode and store the list values
216            tmf_line=tmf_line.replace('!',',')                # get rid of the ! so can use split
217            tmf_line_fields = tmf_line.split(',')             # break the line into fields
218            tmf_line_fields[0] = int(tmf_line_fields[0])      # The tuflow surface RID
219            tmf_line_fields[1] = float(tmf_line_fields[1])    # The tuflow fixed n value
220            tmf_line_fields[2] = int(tmf_line_fields[2])      # Not used
221            tmf_line_fields[3] = int(tmf_line_fields[3])      # Not used
222            tmf_line_fields[4] = float(tmf_line_fields[4])    # The tuflow d1 value
223            tmf_line_fields[5] = float(tmf_line_fields[5])    # The tuflow n1 value
224            tmf_line_fields[6] = float(tmf_line_fields[6])    # The tuflow d2 value
225            tmf_line_fields[7] = float(tmf_line_fields[7])    # The tuflow n2 value
226            tmf_line_fields[8] = str(tmf_line_fields[8])      # The tuflow surface description
227                                       
228            if verbose :
229                print "         ++++ get_tuflow_tmf  appending  "  ,tmf_line_fields
230
231            mat_roughness_data_list.append(tmf_line_fields)   # there are nine data fields including the description
232                       
233   
234               
235    return mat_roughness_data_list            # a list of surface data fields as in the tuflow tmf file
236
237#END FUNCTION ================================================ get_tuflow_tmf
238
239
240#START FUNCTION =============================================== get_tuflow_mat
241def get_tuflow_mat(mat_filename,verbose):
242    'reads in the roughness polys from the tuflow 2d_mat.mif file and IDs from the mid'
243
244    # Read in the polygons from the 2d_mat.mif file
245    matmif_fid=open(mat_filename, 'r')       # open the tuflow material (roughness polys) file
246    matmif_lines=matmif_fid.readlines()      # read all file into string list
247    matmif_fid.close()                       # close the 2d_mat.mif file
248
249   
250    poly_coords_list=[]                               # float list of roughness polygon coord lists
251    point_coords=[]                                   # float list of single coord x y pair
252   
253    for i in range (len(matmif_lines)):
254       
255        if matmif_lines[i][:12].find('Region')>=0: 
256            nlines=int(matmif_lines[i+1] )           # read next line get no of coords (nlines) in poly
257            poly=[]                                  # float list of coords x,y
258           
259            for j in range (nlines):                 # read next nlines of xy coords into poly[]
260                x,y = matmif_lines[i+2+j].split()    # split into two strings
261                x = float(x)                         # convert to numeric coords
262                y = float(y)
263                point_coords = [x,y]
264                poly.append(point_coords)            # add to poly as float xy pairs
265
266            poly_coords_list.append(poly)            # append this poly coord list to list of polys
267            i=i+1+j                                  # reset counter to previous to last line read to avoid problem when at end of file
268
269        else:
270            # unwanted - skip line
271            pass
272     
273
274    # Read in the polygon material (roughness) IDs from the 2d_mat.mid file
275    matmid_fid=open(mat_filename[:-4] + '.mid', 'r') # open the tuflow material (roughness polys) file
276    matmid_lines=matmid_fid.readlines()              # read all file into string list
277    matmid_fid.close()                               # close the 2d_mat.mid file
278
279    poly_RID_list=[]                                 # integer list identifying roughness polygons ID
280    for matmid_line in matmid_lines:
281        RID=int(matmid_line[:(matmid_line.index(','))])  # get the RID of each poly
282        poly_RID_list.append(RID)
283
284    if len(poly_RID_list)== len(poly_coords_list):   # all as it should be
285        if verbose :
286            print "         ++++ get_tuflow_mat - Read in %i surface roughness polygons from %s " %(len(poly_RID_list), mat_filename)
287        return poly_RID_list,poly_coords_list        # [intIDs],[list of poly floatcoords list]
288    else:
289        print "         >>>> ERROR get_tuflow_mat - Poly regions list and poly ID list not of same length in get_tuflow_mat??"
290        print "         >>>> ERROR get_tuflow_mat - Poly ID list %d and poly list %d" %(len(mat_ID_list),len(mat_poly_list))
291        sys.exit(0)                                  # terminate the program
292
293#END FUNCTION ================================================ get_tuflow_mat
294
295
296
Note: See TracBrowser for help on using the repository browser.