""" -----------------------------------------------------------------model_data.py ANUGA FLOOD MODELLING TEMPLATE (Data Assembly) This initialisation script contains surface and scenario-event specific definitions and data for use in the subsequent mesh-domain build process described in model_run.py. This script is imported by the model_run.py script <<<>>> <<<< IT READS IN DATA FROM THE VARIOUS TUFLOW MID/MIF INPUT FILES >>>> Version 1.00: June 2008 Initial release Author : E Rigby (0437 250 500) ted.rigby@rienco.com.au -------------------------------------------------------------------------------- """ ##################################################################################### # # # IMPORT THE REQUIRED PYTHON MODULES # # # ##################################################################################### # Standard python modules import os import time import sys from os import sep # Related ANUGA modules from anuga.utilities.polygon import read_polygon, plot_polygons, polygon_area, is_outside_polygon from anuga.geospatial_data import Geospatial_data from anuga.utilities.system_tools import get_revision_number # Model specific imports import get_tuflow_data # Provides the routines to read in tuflow specific input data ##################################################################################### # # # SET THE REQUIRED INI VARIABLES # # # ##################################################################################### anuga_verbose = True # Ini variable may be set to True to enable additional anuga internal output model_verbose = True # Ini variable may be set to True to enable additional model output model_graphics = False # Ini variable may be set to True to enable model check graphics org_name = 'Rienco Consulting' org_add1 = '4/163 Keira Street' org_add2 = 'Wollongong' org_add3 = 'NSW 2500' org_add4 = 'Australia' org_phone = '0242 250 500' org_web = 'www.rienco.com.au' user_name = 'E.Rigby' user_email = 'ted.rigby@rienco.com.au' ##################################################################################### # # # ASSIGN THE REQUIRED SIMULATION FILES & VARIABLES # # # ##################################################################################### #------------------------------------------------------------------------------------ # Defines the current catchment-simclass-scenario-event and checks files/directories # This assumes a directory structure in which the 'model*'.py files are in a directory # called simclass/scenario. the tms and output directories are checked or created by this # code in a directory below the 'model*'.py files . #------------------------------------------------------------------------------------- # If a historic event simulation (typically for calibration or validation purposes) catchment= 'Macquarie_Rivulet' # Set catchment name simclass = 'Historic' # Set simulation class 'Historic' (actual recorded event) scenario = 'Existing' # Set scenario (eg existing conditions 1991) event = 'Jun91' # Set event (eg Event of June 1991) # If a design event simulation #catchment= 'Macquarie_Rivulet' # Set catchment name #simclass = 'Design' # Set simulation class 'Design' (AR&R design event) #scenario = 'ExAllClear' # Set scenario (eg 'Exallclear','Exblocked','2050AllClear' etc) #event = '100Yr' # Set event ARI (eg 2-100YrARI or PMF) and optionally Duration # Assign the required input files and tms output directory # ANUGA input files # EHR - expand when finalised respoly GIS formats!!!!!!!!!!!! elev_filename = 'macriv_merge_2mdem_trimmed.csv' # assign the raw ALS (csv) filename # ANUGA output/results files results_dir = event # assign output to a directory called 'event' # WBNM - TUFLOW input Files wbnm_filename = 'Macriv_Jun91.wbn' # assign wbnm runfile filename bcdbase_filename = 'bc dbase.csv' # assign tuflow bc datanase filename tmf_filename = 'macriv91.tmf' # assign tuflow tmf filename mat_filename = '2d_mat_Macriv_Jun91.mif' # assign tuflow 2d_mat filename po_filename = '2d_po_Macriv_Jun91.mif' # assign tuflow 2d_po filename # Converted ANUGA tms (hydrographs and dsbc files) directory tms_dir = 'tms_files' # assign anuga tms file directory name # assign a basename to identify screen and result output basename='macriv91' # prefix for output files and screen comments from the simulation # Assign initial water level for this simulation InitialWaterLevel = 0.3 # Set initial lake level to 0.3m AHD # Assign the default RID to apply to the domain (before application of the spatially variable RIDs) default_RID = 2 # This simulation - 'Treed areas of equal resistance at ground and depth' ##################################################################################### # # # CHECK THE SPECIFIED SIMULATION FILES ARE PRESENT # # # ##################################################################################### # Check the required directories and files are present if os.access(tms_dir, os.F_OK) == 0: # Make sure we have the required tms directory below us os.mkdir(tms_dir) print basename+ ' >>>> NOTE - tms directory %s was missing and has been created ' %(tms_dir) if os.access(results_dir, os.F_OK) == 0: # Make sure we have the required results directory below us os.mkdir(results_dir) print basename+ ' >>>> NOTE - Results directory %s was missing and has been created ' %(results_dir) if os.access(elev_filename, os.F_OK) == 0: # check and istop run if absent print basename+ ' >>>> ERROR - specified elevation (ALS) file %s missing ' %(elev_filename) sys.exit(0) # terminate the program if os.access(wbnm_filename, os.F_OK) == 0: # check and istop run if absent print basename+ ' >>>> ERROR - specified wbnm runfile %s missing ' %(wbnm_filename) sys.exit(0) # terminate the program if os.access(bcdbase_filename, os.F_OK) == 0: # check and stop run if absent print basename+ ' >>>> ERROR - specified tuflow bcdbase file %s missing ' %(bcdbase_filename) sys.exit(0) # terminate the program if os.access(tmf_filename, os.F_OK) == 0: # check and stop run if absent print basename+ ' >>>> ERROR - specified tuflow btmf file %s missing ' %(tmf_filename) sys.exit(0) # terminate the program if os.access(mat_filename, os.F_OK) == 0: # check and stop run if absent print basename+ ' >>>> ERROR - specified tuflow 2d_mat file %s missing ' %(mat_filename) sys.exit(0) # terminate the program if os.access(po_filename, os.F_OK) == 0: # check and stop run if absent print basename+ ' >>>> ERROR - specified tuflow 2d_po file %s missing ' %(po_filename) sys.exit(0) # terminate the program print basename+ ' >>>> #########################################################################' print basename+ ' >>>> COMMENCING THIS MODEL SIMULATION USING ANUGA REVISION NUMBER %s ' %0#%(get_revision_number()) print basename+ ' >>>> for the simulation named ' print basename+ ' >>>> ',catchment,'-',simclass,'-',scenario,'-',event print basename+ ' >>>> Run dated %s by %s of %s ' %(time.ctime(),user_name,org_name) print basename+ ' >>>> #########################################################################\n' print basename+ ' >>>> Checked - All required files/directories exist - OK to proceed! \n ' ##################################################################################### # # # COMMENCE ASSEMBLY OF THE REQUIRED DATA # # # ##################################################################################### print basename+' >>>> Commencing data assembly for this simulation ' ##################################################################################### # # # ASSEMBLE THE BOUNDING POLYGON # # # ##################################################################################### #------------------------------------------------------------------------------ # Define bounding polygon and default res to be used later in model_run.py to build # the mesh and boudary strings with their tags # # OUTPUTS OF THIS CODE BLOCK ARE; # bounding_polygon -- The bounding polygon within which the mesh will be created # default_res -- The default resolution to apply in creating the mesh #------------------------------------------------------------------------------ print basename+ ' >>>> Assemble the bounding polygon and set default mesh resolution' # Note..Bounding poly is a single contiguous polygon initially created in MapInfo. # Once selected it was exported to MID/MIF and a bounding_poly,csv created by editing the MIF # to remove pre and postambles and to add comma separators # EHR...Clumsy - need to be able to read bounding poly directly from (tuflow?)GIS!!! bounding_polygon = read_polygon('bounding_poly.csv') print basename+ ' >>>> Area of bounding polygon in Km2', polygon_area(bounding_polygon)/1000000.0 # Set mesh filename to be created with *.py files in simclass/scenario directory mesh_filename = basename+'.msh' # Set bounding polygon (default) mesh resolution default_res = 200 # m^2 (max area of triangle in triangular mesh) (coarse res rural and forest)) print basename+ ' >>>> Default resolution within bounding polygon is ',default_res,' m2' ##################################################################################### # # # ASSEMBLE THE TAGGED BOUNDARY STRINGS AROUND THE BOUNDING POLYGON # # # ##################################################################################### #------------------------------------------------------------------------------------ # Create tagged boundary string list for use in create-mesh in model_build.py # EHR... finding out point IDs to use is a hassle . Used mapinfo to identify vertex # numbers and then manually transfered to model_data.py (dont forget to subtract 1??) # Needs some way of creating in GIS and being read in?? # # OUTPUTS OF THIS CODE BLOCK ARE; # boundary_strings[]-- List of vertices of each tagged string segment around the # bounding polygon at which boundary conditions will be defined #------------------------------------------------------------------------------------ print basename+ ' >>>> Assemble tagged boundary strings around the bounding polygon ' boundary_strings = {'Lake': [84,85,86], 'MaMt': [56], 'Mriv': [39], 'Cbck': [36], 'WFra': [20], 'EFra': [12,13] } # Note Anuga will automatically create boundary tagged strings for what is left of the bounding # poly with the tag 'exterior'. ALL must be associated with a BC in the mainline Model_run.py print basename+ ' >>>> %i tagged boundary strings assigned by user' %(len(boundary_strings)) ##################################################################################### # # # ASSEMBLE INTERIOR RESOLUTION POLYGONS # # # ##################################################################################### #------------------------------------------------------------------------------- # Define Interior region polygons for later use in model_run.py to define zones # of different mesh resolution # This is an ANUGA specific input with no Tuflow equivalent # # OUTPUTS OF THIS CODE BLOCK ARE; # interior_resregions[]-- List of interior res regions containing [res_poly, res] # to be used in model_run to create the mesh #------------------------------------------------------------------------------- print basename+ ' >>>> Assemble any inner region polygons and set associated resolutions' # EHR...Clumsy - need to be able to read multiple regions directly from GIS fine_polygon = read_polygon('fine_poly.csv') # Fine) res (typ streams and around bridge openingss) fine_res = 50 # m^2 max #Note coarse not needed as set as default for bounding_polgon # Plot res polygons out to event/plot_respoly_filename.png for visual check if model_graphics : plot_respoly_filename = event+sep+'Macriv_Res_Polys' plot_polygons([bounding_polygon, fine_polygon], style=None,figname=plot_respoly_filename,verbose=model_verbose) # Create interior region resolutions list for later use in model_run.py - only one for this model interior_resregions = [ [ fine_polygon, fine_res] ] print basename+ ' >>>> %i interior resolution regions added ' %(len(interior_resregions)) ##################################################################################### # # # ASSEMBLE ELEVATION DATA # # # ##################################################################################### #------------------------------------------------------------------------------ # Assign reference to file with elevation data for later uset in model_build.py # the following assumes these elevation data files are stored with the *.py files # in the simclass/scenario (current) directory # # OUTPUTS OF THIS CODE BLOCK ARE; # elev_filename -- Filename containing the ALS data to be used to set mesh elevations #------------------------------------------------------------------------------ print basename+ ' >>>> Assemble the file containing elevation (ALS) data from %s ' %(elev_filename) # Assign/pre-process the file with the data to be used to assign mesh elevations # Identify the raw ALS data file # OPTION 1 (Preferred but limited to smaller ALS datasets) # This option converts the raw xyz ALS data into NETCDF (pts) format # This only needs to be done once. The pts file is then used to assign # mesh elevations in model_run.py # #print basename+' >>>> Elevation x,y,z scatterpoint dataset %s being used to create NETCDF pts file' %(elev_filename) # Convert into NETCDF pts format to speed later use in setting mesh elevations #G = Geospatial_data(elev_filename, max_read_lines=500, load_file_now=True, verbose=True) #G.export_points_file(basename+'.pts') # switch elev_filename to now point to the NETCDF pts file #elev_filename = basename+'.pts' # OPTION 2 (Leave elevation data in original (csv) format) # Export_points loads all data into memory (no blocking) so is limited in # regard to the number of points it can read in/export. # If there are more than about 6 million points in the ALS dataset # export_points will fail with a memory over-run error. # In such a case the elevation data should be left in csv format. # The elevation data is left in csv format in this model as the ALS dataset overruns memory ! print basename+ ' >>>> %s will be used to supply elevation data to the mesh' %(elev_filename) ##################################################################################### # # # ASSEMBLE TUFLOW SURFACE ROUGHNESS POLYGONS AND ASSOCIATED DATA # # # ##################################################################################### #------------------------------------------------------------------------------- # Read in tuflow roughness polygons and associated roughness data for later use # in model_run.py # # OUTPUTS OF THIS CODE BLOCK ARE; # mat_roughness_data[RID] -- dictionary of properties for each surface (RID) # mat_RID_list[RID] -- surface roughness Ids (RID) for each poly # mat_poly_list[[x0,y0],[x1,y1]...] -- list of vertices for each roughness poly #-------------------------------------------------------------------------------- print basename+ ' >>>> Assemble the surface roughness data from Tuflows %s and %s ' %(tmf_filename, mat_filename) # Get the tuflow roughness data for each surface class mat_roughness_data_list=[] mat_roughness_data_list=get_tuflow_data.get_tuflow_tmf(tmf_filename,verbose=model_verbose) # Turn material roughness data list into a dictionary mat_roughness_data = {} for entry in mat_roughness_data_list: key = entry[0] # rid mat_roughness_data[key] = [entry[1], # n0 entry[4], # d1 entry[5], # n1 entry[6], # d2 entry[7]] # n2 mat_RID_list=[] # Int list of rough poly IDs-one for each poly mat_poly_list=[] # list of float list of poly coords mat_RID_list,mat_poly_list=get_tuflow_data.get_tuflow_mat(mat_filename,model_verbose) # get the tuflow roughness polys # Plot roughness polygons out to event/plot_roughpoly_filename.png for visual check if model_graphics : plot_roughpoly_filename = event+sep+'Macriv_rough_Polys' plot_polygons( mat_poly_list, style=None, figname=plot_roughpoly_filename,verbose=model_verbose) # Extract the default RID roughness data for i in range(len(mat_roughness_data_list)) : if default_RID == mat_roughness_data_list[i][0] : # got a match default_n = mat_roughness_data_list[i][1] # default_n is initially the tuflow fixed n value default_mat = mat_roughness_data_list[i][8].rstrip(' ') # default mat name if model_verbose: print ' ++++ The default RID is %i - fixed n=%3.3f - %s ' %(default_RID,default_n,default_mat) print basename+ ' >>>> %i roughness polys and %i unique roughness IDs read in' %(len(mat_poly_list), len(mat_roughness_data_list)) ####################################################################################### # # # ASSEMBLE TUFLOW (TEMPORAL) SUBAREAL INFLOWS AND DOWNSTREAM BOUNDARY CONDITIONS # # # ####################################################################################### #-------------------------------------------------------------------------------------- # This section reads in the tuflow bc_dbase.csv file to obtain locations and file # references for the associated inflows or downstream boundary conditions. # It then creates the NetCDF flow hydrograph files from the tuflow ts1 file. # model_run then uses these to load the local or total temporal inflows into the domain. # The downstream boundary conditions are read in at the same time and stored for # later application in model_run where boundary conditions are applied. # # OUTPUTS OF THIS CODE BLOCK ARE; # inflow_hydrographs[]-- List of inflow hydrographs containing [tms_filename, xcoord,ycoord] # dsbc_hydrographs[] -- List of dsbc hydrographs containing [dsbc_filename] #-------------------------------------------------------------------------------------- # ASSEMBLE THE (TEMPORAL) SUBAREAL INFLOWS::::::::::::::::::::::::::::::::::::::::::::::: # Note: The current inflow function 'pours' flow onto the domain within a circle or # polygonal area specified by the user at the specified temporal rate. # In due course, the above inflows should be augmented by latteral inflows at the boundary # (total as distinct from local hydrographs) to better correspond with conventional # flood modelling. # A forcing function for flow (or rainfall) poured onto a polygon is still however # relevant to inflow/rainfall falling directly onto the model domain print basename+" >>>> Assemble the inflow and dsbc hydrographs from the tuflow database file", bcdbase_filename inflowNo = int(-1) # Initialise no of inflow hydrographs dsbcNo = int(-1) # Initialise no of dsbc hydrographs inflow_hydrographs = [] # Initialise list of inflow hydros dsbc_hydrograph = [] # Initialise dsbc hydrograph inflow_hydrograph = [] # Initialise inflow hydrograph outlet = [] # initialise x,y coords of subareal outlet dsbc_hydrographs = [] # Initialise list of dsbc hydrographs bc_fid=open(bcdbase_filename, 'r') # open the primary bc_dbase file bc_lines = bc_fid.readlines() # read all file into string list bc_fid.close() # close the bc_dbase file for bc_line in bc_lines[1:]: # skip headers and read each line of bc_lines and process # break each bcdbase line into fields confirm type of bc data and obtain/append hydrograph bc_fields=bc_line.split(',') # split bc_line into bc fields location=bc_fields[0] # record the bc location loc_tempname=bc_fields[1] # record the assoc bc's filename loc_filename=loc_tempname.replace("_EVENT_",'_'+event+'_') # modify loc filename to reflect event bc_filetype=get_tuflow_data.get_tuflow_bc_filetype(loc_filename,verbose=model_verbose) # obtain type = "ts1"(flowrate) or "csv" (dsbc elev) if loc_filename.find(".ts1")>=0: # Is file containing hydrographs for each subarea bc_filetype="ts1" elif loc_filename.find(".csv")>=0: # Is time(hrs),elev(m) dsbc.csv file bc_filetype="csv" else: print basename+ " >>>> ERROR model_data.py - Unrecognised bc filetype %s - aborting " %(bc_filetype) sys.exit(0) # terminate the program if bc_filetype == "ts1": # add inflow hydrograph inflowNo = inflowNo +1 # increment inflow hydro Id # read in the inflow time series hyd_time,hyd_flow=get_tuflow_data.get_tuflow_bc(location=location, loc_filename=loc_filename, bc_filetype=bc_filetype, verbose=model_verbose) # convert to tms NetCDF file type tms_filename=get_tuflow_data.convert_asc_ts_to_tms(tms_dir,location,'seconds',hyd_time,hyd_flow,verbose=model_verbose) # get coords of subareal outlet node xcoord,ycoord=get_tuflow_data.get_wbnm_coords(wbnm_filename=wbnm_filename,subarea=location,verbose=model_verbose) outlet=[xcoord,ycoord] # check inflow coords are within bounding poly if is_outside_polygon(outlet,bounding_polygon,closed = True,verbose=False): 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) sys,exit(0) # terminate the program inflow_hydrograph=[tms_filename,xcoord,ycoord] inflow_hydrographs.append(inflow_hydrograph) # Store details ready to use in model_run if bc_filetype == "csv": # add dsbc hydrograph dsbcNo = dsbcNo+1 # Increment dsbc hydro ID # read in downstream boundary condition water elevation time series dsbc_time,dsbc_elev=get_tuflow_data.get_tuflow_bc(location=location, loc_filename=loc_filename, bc_filetype=bc_filetype,verbose=model_verbose) # convert to tms NetCDF file type dsbc_tms_filename=get_tuflow_data.convert_asc_ts_to_tms(tms_dir,location,'seconds',dsbc_time,dsbc_elev,verbose=model_verbose) dsbc_hydrographs.append(dsbc_tms_filename) # Store details ready to use in model_run if dsbcNo > 1 : print basename+ " >>>> ERROR model_data_py - Current model is limited to 1 but dsbc - %i read in" %(dsbcNo) sys.exit(0) # EHR -- need to generalise to many dsbc!!!!!! print basename+' >>>> %i inflow and %i dsbc hydrographs read in and and stored ' %(len(inflow_hydrographs),len(dsbc_hydrographs)) print basename+' >>>> Converted tuflow inflow and dsbc hydrographs and stored as tms files in %s ' %(sep+tms_dir)