source: anuga_work/development/ted_rigby/model_data.py @ 5973

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

Ted Rigby's variable friction code

File size: 25.5 KB
Line 
1"""
2-----------------------------------------------------------------model_data.py
3ANUGA FLOOD MODELLING TEMPLATE (Data Assembly)
4This initialisation script contains surface and scenario-event specific definitions and
5data for use in the subsequent mesh-domain build process described in model_run.py.
6This script is imported by the model_run.py script
7
8<<<<THIS VERSION HAS BEEN CONSTRUCTED TO ACCEPT TUFLOW DATA INPUTS>>>>
9<<<< IT READS IN DATA FROM THE VARIOUS TUFLOW MID/MIF INPUT FILES >>>>
10
11Version 1.00: June 2008 Initial  release
12Author      : E Rigby (0437 250 500) ted.rigby@rienco.com.au
13
14--------------------------------------------------------------------------------
15"""
16#####################################################################################
17#                                                                                   #
18#                     IMPORT THE REQUIRED PYTHON MODULES                            #
19#                                                                                   #
20#####################################################################################
21
22# Standard python modules
23import os
24import time
25import sys
26from os import sep
27
28# Related ANUGA modules
29from anuga.utilities.polygon import read_polygon, plot_polygons, polygon_area, is_outside_polygon
30from anuga.geospatial_data import Geospatial_data
31from anuga.utilities.system_tools import get_revision_number
32
33# Model specific imports
34import get_tuflow_data      # Provides the routines to read in tuflow specific input data
35
36#####################################################################################
37#                                                                                   #
38#                       SET THE REQUIRED INI VARIABLES                              #
39#                                                                                   #
40#####################################################################################
41
42anuga_verbose   = True       # Ini variable may be set to True to enable additional anuga internal output
43model_verbose   = True       # Ini variable may be set to True to enable additional model output
44model_graphics  = False      # Ini variable may be set to True to enable model check graphics
45
46org_name        = 'Rienco Consulting'
47org_add1        = '4/163 Keira Street'
48org_add2        = 'Wollongong'
49org_add3        = 'NSW 2500'
50org_add4        = 'Australia'
51org_phone       = '0242 250 500'
52org_web         = 'www.rienco.com.au'
53
54user_name       = 'E.Rigby'
55user_email      = 'ted.rigby@rienco.com.au'
56
57
58#####################################################################################
59#                                                                                   #
60#            ASSIGN  THE REQUIRED SIMULATION FILES & VARIABLES                      #
61#                                                                                   #
62#####################################################################################
63
64#------------------------------------------------------------------------------------
65# Defines the  current catchment-simclass-scenario-event  and checks files/directories
66# This assumes a directory structure in which the  'model*'.py files are in a directory
67# called simclass/scenario. the tms and output directories are checked or created by this
68# code in  a directory below the 'model*'.py files .
69#-------------------------------------------------------------------------------------
70
71
72# If a historic event simulation (typically for calibration or validation purposes)
73catchment= 'Macquarie_Rivulet'               # Set catchment name
74simclass = 'Historic'                        # Set simulation class 'Historic' (actual recorded event)
75scenario = 'Existing'                        # Set scenario (eg existing conditions 1991)
76event    = 'Jun91'                           # Set event (eg Event of June 1991)
77
78# If a design event simulation
79#catchment= 'Macquarie_Rivulet'               # Set catchment name
80#simclass = 'Design'                          # Set simulation class  'Design' (AR&R design event)
81#scenario = 'ExAllClear'                      # Set scenario (eg 'Exallclear','Exblocked','2050AllClear' etc)
82#event    = '100Yr'                           # Set event ARI (eg 2-100YrARI or PMF) and optionally Duration
83
84# Assign the required input files and tms output directory
85
86# ANUGA input files
87# EHR - expand when finalised respoly GIS formats!!!!!!!!!!!!
88elev_filename     = 'macriv_merge_2mdem_trimmed.csv' # assign the raw ALS (csv) filename
89
90# ANUGA output/results files
91results_dir       = event                     # assign output to a directory called 'event'
92
93# WBNM - TUFLOW input Files
94wbnm_filename     = 'Macriv_Jun91.wbn'        # assign wbnm runfile filename
95bcdbase_filename  = 'bc dbase.csv'            # assign tuflow bc datanase filename
96tmf_filename      = 'macriv91.tmf'            # assign tuflow tmf filename
97mat_filename      = '2d_mat_Macriv_Jun91.mif' # assign tuflow 2d_mat filename
98po_filename       = '2d_po_Macriv_Jun91.mif'  # assign tuflow 2d_po filename
99
100# Converted ANUGA tms (hydrographs and dsbc files) directory
101tms_dir           = 'tms_files'               # assign anuga tms file directory name
102
103# assign a basename to identify screen and result output
104basename='macriv91'                           # prefix for output files and screen comments from the simulation
105
106# Assign initial water level for this simulation
107InitialWaterLevel = 0.3                       # Set initial lake level to 0.3m AHD
108
109# Assign the default RID to apply to the domain (before application of the spatially variable RIDs)
110default_RID = 2                               #  This simulation - 'Treed areas of equal resistance at ground and depth'
111
112
113#####################################################################################
114#                                                                                   #
115#            CHECK  THE SPECIFIED SIMULATION FILES ARE PRESENT                      #
116#                                                                                   #
117#####################################################################################
118
119# Check the required  directories and files are present
120if os.access(tms_dir, os.F_OK) == 0:          # Make sure we have the required tms directory below us
121    os.mkdir(tms_dir)
122    print basename+ ' >>>> NOTE - tms directory %s was missing and has been created  ' %(tms_dir)
123   
124if os.access(results_dir, os.F_OK) == 0:      # Make sure we have the required  results directory below us
125    os.mkdir(results_dir)
126    print basename+ ' >>>> NOTE - Results directory %s was missing and has been created  ' %(results_dir)
127
128if os.access(elev_filename, os.F_OK) == 0:    # check and istop run if absent
129    print basename+ ' >>>> ERROR - specified elevation (ALS) file %s missing ' %(elev_filename)
130    sys.exit(0)                               # terminate the program
131   
132if os.access(wbnm_filename, os.F_OK) == 0:    # check and istop run if absent
133    print basename+ ' >>>> ERROR - specified wbnm runfile %s missing ' %(wbnm_filename)
134    sys.exit(0)                               # terminate the program
135   
136if os.access(bcdbase_filename, os.F_OK) == 0: # check and stop run if absent
137    print basename+ ' >>>> ERROR - specified tuflow bcdbase file %s missing ' %(bcdbase_filename)
138    sys.exit(0)                               # terminate the program
139   
140if os.access(tmf_filename, os.F_OK) == 0:     # check and stop run if absent
141    print basename+ ' >>>> ERROR - specified tuflow btmf file %s missing ' %(tmf_filename)
142    sys.exit(0)                               # terminate the program
143   
144if os.access(mat_filename, os.F_OK) == 0:     # check and stop run if absent
145    print basename+ ' >>>> ERROR - specified tuflow 2d_mat file %s missing ' %(mat_filename)
146    sys.exit(0)                               # terminate the program
147   
148if os.access(po_filename, os.F_OK) == 0:      # check and stop run if absent
149    print basename+ ' >>>> ERROR - specified tuflow 2d_po file %s missing ' %(po_filename)
150    sys.exit(0)                               # terminate the program
151
152
153print basename+ ' >>>> #########################################################################'
154print basename+ ' >>>>     COMMENCING THIS MODEL SIMULATION USING ANUGA REVISION NUMBER %s  ' %0#%(get_revision_number())
155print basename+ ' >>>>                         for the simulation named        '
156print basename+ ' >>>>             ',catchment,'-',simclass,'-',scenario,'-',event
157print basename+ ' >>>>     Run dated %s by %s of %s ' %(time.ctime(),user_name,org_name)
158print basename+ ' >>>> #########################################################################\n'
159
160print basename+ ' >>>> Checked - All required files/directories exist - OK to proceed! \n '
161
162#####################################################################################
163#                                                                                   #
164#                    COMMENCE ASSEMBLY OF THE REQUIRED DATA                         #
165#                                                                                   #
166#####################################################################################
167   
168print basename+' >>>> Commencing data assembly for this simulation '
169
170
171#####################################################################################
172#                                                                                   #
173#                        ASSEMBLE THE BOUNDING POLYGON                              #
174#                                                                                   #
175#####################################################################################
176
177#------------------------------------------------------------------------------
178# Define bounding polygon and default res to be used later in model_run.py to build
179# the mesh and boudary strings with their tags
180#
181# OUTPUTS OF THIS CODE BLOCK ARE;
182# bounding_polygon -- The bounding polygon within which the mesh will be created
183# default_res      -- The default resolution to apply in creating the mesh
184#------------------------------------------------------------------------------
185
186print basename+ ' >>>> Assemble the bounding polygon and set default mesh resolution'
187
188
189# Note..Bounding poly is a single contiguous polygon initially created in MapInfo.
190#       Once selected it was exported to MID/MIF and a bounding_poly,csv created by editing the MIF
191#       to remove pre and postambles and to add comma separators
192# EHR...Clumsy - need to be able to read bounding poly directly from (tuflow?)GIS!!!
193bounding_polygon = read_polygon('bounding_poly.csv')
194
195print basename+ ' >>>> Area of bounding polygon in Km2', polygon_area(bounding_polygon)/1000000.0
196
197# Set mesh filename to be created with  *.py files in simclass/scenario directory
198mesh_filename = basename+'.msh'
199
200# Set bounding polygon (default)  mesh resolution
201default_res = 200      # m^2  (max area of triangle in  triangular mesh) (coarse res rural and forest))
202
203print basename+ ' >>>> Default resolution within bounding polygon is ',default_res,' m2'
204
205
206#####################################################################################
207#                                                                                   #
208#         ASSEMBLE THE TAGGED BOUNDARY STRINGS AROUND THE BOUNDING POLYGON          #
209#                                                                                   #
210#####################################################################################
211
212#------------------------------------------------------------------------------------
213# Create  tagged boundary string list for use in create-mesh in model_build.py
214# EHR... finding out point IDs to use is a hassle . Used mapinfo to identify vertex
215#        numbers and then manually transfered to model_data.py (dont forget to subtract 1??)
216#        Needs some way of creating in GIS and being read in??
217#
218# OUTPUTS OF THIS CODE BLOCK ARE;
219# boundary_strings[]-- List of vertices of each tagged string segment around the 
220#                      bounding polygon at which boundary conditions will be defined
221#------------------------------------------------------------------------------------
222
223print basename+ ' >>>> Assemble tagged boundary strings around the bounding polygon '
224
225boundary_strings = {'Lake': [84,85,86],
226                    'MaMt': [56],
227                    'Mriv': [39],
228                    'Cbck': [36],
229                    'WFra': [20],
230                    'EFra': [12,13]  }
231
232# Note Anuga will automatically create boundary tagged strings for what is left of the bounding
233# poly with the tag 'exterior'. ALL must be associated with a BC in the mainline Model_run.py
234
235print basename+ ' >>>> %i tagged boundary strings assigned by user' %(len(boundary_strings))
236
237#####################################################################################
238#                                                                                   #
239#                   ASSEMBLE INTERIOR RESOLUTION POLYGONS                           #
240#                                                                                   #
241#####################################################################################
242
243#-------------------------------------------------------------------------------
244# Define Interior region polygons for later use in model_run.py to define zones
245# of different mesh resolution
246# This is an ANUGA specific input with no Tuflow equivalent
247#
248# OUTPUTS OF THIS CODE BLOCK ARE;
249# interior_resregions[]-- List of interior res regions containing [res_poly, res]
250#                         to be used in model_run to create the mesh
251#-------------------------------------------------------------------------------
252
253print basename+ ' >>>> Assemble any inner region polygons and set associated resolutions'
254
255# EHR...Clumsy - need to be able to read multiple regions directly from GIS
256fine_polygon      = read_polygon('fine_poly.csv')              # Fine) res (typ streams and around bridge openingss)
257fine_res = 50                                                  # m^2 max
258#Note coarse not needed as  set as default for bounding_polgon
259
260# Plot res polygons out to event/plot_respoly_filename.png for visual check
261if model_graphics :
262    plot_respoly_filename = event+sep+'Macriv_Res_Polys'
263    plot_polygons([bounding_polygon, fine_polygon], style=None,figname=plot_respoly_filename,verbose=model_verbose)
264
265# Create interior region resolutions list for later use in model_run.py - only one for this model
266interior_resregions = [ [ fine_polygon, fine_res] ]
267
268print basename+ ' >>>> %i interior resolution regions added ' %(len(interior_resregions))
269
270#####################################################################################
271#                                                                                   #
272#                           ASSEMBLE ELEVATION DATA                                 #
273#                                                                                   #
274#####################################################################################
275
276#------------------------------------------------------------------------------
277# Assign reference to file with elevation data for later uset in model_build.py
278# the following assumes these elevation data files are stored with the *.py files
279# in the simclass/scenario (current) directory
280#
281# OUTPUTS OF THIS CODE BLOCK ARE;
282# elev_filename -- Filename containing the ALS data to be used to set mesh elevations
283#------------------------------------------------------------------------------
284
285print basename+ ' >>>> Assemble the file containing elevation (ALS) data from %s ' %(elev_filename)
286
287# Assign/pre-process the file with the data to be used to assign mesh elevations
288# Identify the raw ALS data file
289
290# OPTION 1 (Preferred but limited to smaller ALS datasets)
291#           This option converts the raw xyz ALS data into NETCDF (pts) format
292#           This only needs to be done once. The pts file is then used to assign
293#           mesh elevations in model_run.py
294#
295#print basename+' >>>> Elevation x,y,z scatterpoint dataset %s being used to create NETCDF pts file' %(elev_filename)
296# Convert into NETCDF pts format to speed later use in setting mesh elevations
297#G = Geospatial_data(elev_filename, max_read_lines=500, load_file_now=True, verbose=True)
298#G.export_points_file(basename+'.pts')
299# switch elev_filename to now point to the NETCDF pts file
300#elev_filename = basename+'.pts'
301
302# OPTION 2 (Leave elevation data in original (csv) format)
303#           Export_points loads all data into memory (no blocking) so is limited in
304#           regard to the number of points it can read in/export.
305#           If there are more than about 6 million points in the ALS dataset
306#           export_points will fail with a memory over-run error.
307#           In such a case the elevation data should be left in csv format.
308
309# The elevation data is left in csv format in this model as the ALS dataset overruns memory !
310
311print basename+ ' >>>> %s will be used to supply elevation data to the mesh' %(elev_filename)
312
313#####################################################################################
314#                                                                                   #
315#         ASSEMBLE TUFLOW SURFACE ROUGHNESS POLYGONS AND ASSOCIATED DATA            #
316#                                                                                   #
317#####################################################################################
318
319#-------------------------------------------------------------------------------
320# Read in tuflow roughness polygons and associated roughness data for later use
321# in model_run.py
322#
323# OUTPUTS OF THIS CODE BLOCK ARE;
324# mat_roughness_data_list[[RID],[n0] etc.] -- list of properties for each surface (RID)
325# mat_RID_list[RID]                        --  surface roughness Ids (RID) for each poly
326# mat_poly_list[[x0,y0],[x1,y1]...]        --  list of vertices for each roughness poly
327#--------------------------------------------------------------------------------
328
329
330print basename+ ' >>>> Assemble the surface roughness data from Tuflows %s and %s ' %(tmf_filename, mat_filename)
331
332mat_roughness_data_list=[]
333mat_roughness_data_list=get_tuflow_data.get_tuflow_tmf(tmf_filename,verbose=model_verbose)    # get the tuflow roughness data for each surface class
334
335mat_RID_list=[]                                                         # Int list of rough poly IDs-one for each poly
336mat_poly_list=[]                                                        # list of float list of poly coords
337
338mat_RID_list,mat_poly_list=get_tuflow_data.get_tuflow_mat(mat_filename,model_verbose) # get the tuflow roughness polys
339
340# Plot  roughness polygons out to event/plot_roughpoly_filename.png for visual check
341if model_graphics :
342    plot_roughpoly_filename = event+sep+'Macriv_rough_Polys'
343    plot_polygons( mat_poly_list, style=None,  figname=plot_roughpoly_filename,verbose=model_verbose) 
344
345# Extract the default RID roughness data
346for i in range(len(mat_roughness_data_list)) :
347    if default_RID == mat_roughness_data_list[i][0] :           # got a match
348        default_n = mat_roughness_data_list[i][1]               # default_n is initially the tuflow fixed n value
349        default_mat = mat_roughness_data_list[i][8].rstrip(' ') # default mat name
350       
351if model_verbose:
352    print '         ++++ The default RID is %i - fixed n=%3.3f - %s ' %(default_RID,default_n,default_mat)
353
354print basename+ ' >>>> %i roughness polys and %i unique roughness IDs read in' %(len(mat_poly_list), len(mat_roughness_data_list))
355
356
357#######################################################################################
358#                                                                                     #
359#  ASSEMBLE TUFLOW (TEMPORAL) SUBAREAL INFLOWS AND DOWNSTREAM BOUNDARY CONDITIONS     #
360#                                                                                     #
361#######################################################################################
362
363#--------------------------------------------------------------------------------------
364# This section reads in the tuflow bc_dbase.csv file to obtain locations and file
365# references for the associated inflows or downstream boundary conditions.
366# It then creates the NetCDF flow hydrograph files from the tuflow ts1 file.
367# model_run then uses these to load the local or total temporal inflows into the domain.
368# The downstream boundary conditions are read in at the same time and stored for
369# later application in model_run where boundary conditions are applied.
370#
371# OUTPUTS OF THIS CODE BLOCK ARE;
372# inflow_hydrographs[]-- List of inflow hydrographs containing [tms_filename, xcoord,ycoord]
373# dsbc_hydrographs[]  -- List of   dsbc hydrographs containing [dsbc_filename]
374#--------------------------------------------------------------------------------------
375
376
377# ASSEMBLE THE (TEMPORAL) SUBAREAL INFLOWS:::::::::::::::::::::::::::::::::::::::::::::::
378
379# Note: The current inflow function 'pours' flow onto the domain within a circle or
380# polygonal area specified by the user at the specified temporal rate.
381# In due course, the above inflows should be augmented by latteral inflows at the boundary
382# (total as distinct from local hydrographs) to better correspond with conventional
383# flood modelling.
384# A forcing function for flow (or rainfall) poured onto a polygon is still however
385# relevant to inflow/rainfall falling directly onto the model domain
386
387print basename+" >>>> Assemble the inflow and dsbc hydrographs from the tuflow database file", bcdbase_filename
388
389inflowNo = int(-1)                            # Initialise no of inflow hydrographs
390dsbcNo   = int(-1)                            # Initialise no of dsbc hydrographs
391inflow_hydrographs = []                       # Initialise list of inflow hydros
392dsbc_hydrograph    = []                       # Initialise dsbc hydrograph
393inflow_hydrograph  = []                       # Initialise inflow hydrograph
394outlet             = []                       # initialise x,y coords of subareal outlet
395dsbc_hydrographs   = []                       # Initialise list of dsbc hydrographs
396
397bc_fid=open(bcdbase_filename, 'r')            # open the primary bc_dbase file
398bc_lines = bc_fid.readlines()                 # read all file into string list
399bc_fid.close()                                # close the bc_dbase file
400     
401for bc_line in bc_lines[1:]:                  # skip headers and read each line of bc_lines and process
402   
403    # break each bcdbase line into fields  confirm type of bc data and obtain/append hydrograph
404   
405    bc_fields=bc_line.split(',')             # split bc_line into bc fields
406    location=bc_fields[0]                    # record the bc location
407    loc_tempname=bc_fields[1]                # record the assoc bc's filename
408    loc_filename=loc_tempname.replace("_EVENT_",'_'+event+'_') # modify loc filename to reflect event
409   
410    bc_filetype=get_tuflow_data.get_tuflow_bc_filetype(loc_filename,verbose=model_verbose)   # obtain type = "ts1"(flowrate) or "csv" (dsbc elev)
411
412    if loc_filename.find(".ts1")>=0:         # Is file containing hydrographs for each subarea
413       bc_filetype="ts1"
414               
415    elif loc_filename.find(".csv")>=0:       # Is time(hrs),elev(m) dsbc.csv file
416        bc_filetype="csv"
417               
418    else:
419        print basename+ " >>>> ERROR model_data.py - Unrecognised bc filetype %s - aborting " %(bc_filetype)
420        sys.exit(0)                          # terminate the program
421 
422    if bc_filetype == "ts1":                 # add inflow hydrograph
423        inflowNo = inflowNo +1               # increment inflow hydro Id
424        # read in the inflow time series
425        hyd_time,hyd_flow=get_tuflow_data.get_tuflow_bc(location=location, loc_filename=loc_filename, bc_filetype=bc_filetype, verbose=model_verbose)   
426        # convert to tms NetCDF file type
427        tms_filename=get_tuflow_data.convert_asc_ts_to_tms(tms_dir,location,'seconds',hyd_time,hyd_flow,verbose=model_verbose)
428        # get coords of subareal outlet node
429        xcoord,ycoord=get_tuflow_data.get_wbnm_coords(wbnm_filename=wbnm_filename,subarea=location,verbose=model_verbose)
430        outlet=[xcoord,ycoord]
431        # check inflow coords are within bounding poly
432        if is_outside_polygon(outlet,bounding_polygon,closed = True,verbose=False):
433            print basename+ " >>>> ERROR model_data.py - Subareal outlet for %s at - %7.2f %7.2f is on or outside the bounding polygon " %(tms_filename,xcoord,ycoord)
434            sys,exit(0)                     # terminate the program
435                       
436        inflow_hydrograph=[tms_filename,xcoord,ycoord]
437        inflow_hydrographs.append(inflow_hydrograph)   # Store details ready to use in model_run
438       
439       
440    if bc_filetype == "csv":                # add dsbc hydrograph
441        dsbcNo = dsbcNo+1                   # Increment dsbc hydro ID
442        # read in downstream boundary condition water elevation time series
443        dsbc_time,dsbc_elev=get_tuflow_data.get_tuflow_bc(location=location, loc_filename=loc_filename, bc_filetype=bc_filetype,verbose=model_verbose)   
444        # convert to tms NetCDF file type
445        dsbc_tms_filename=get_tuflow_data.convert_asc_ts_to_tms(tms_dir,location,'seconds',dsbc_time,dsbc_elev,verbose=model_verbose)
446        dsbc_hydrographs.append(dsbc_tms_filename)    # Store details ready to use in model_run           
447        if dsbcNo > 1 :
448            print basename+ " >>>> ERROR model_data_py - Current model is limited to 1 but dsbc - %i read in" %(dsbcNo)
449            sys.exit(0)                               # EHR -- need to generalise to many dsbc!!!!!!
450                       
451print basename+' >>>> %i inflow  and %i dsbc hydrographs read in and and stored ' %(len(inflow_hydrographs),len(dsbc_hydrographs)) 
452
453print basename+' >>>> Converted tuflow inflow  and dsbc hydrographs and stored as tms files in %s ' %(sep+tms_dir)
454
455
Note: See TracBrowser for help on using the repository browser.