""" -----------------------------------------------------------------model_data.py ANUGA TSUNAMI MODELLING TEMPLATE (Data Assembly) This script extracts/converts client model-scenario-event data for use in the subsequent Anuga mesh-domain build process described in model_build.py. Version 1.00 June 2008 Initial release for tuflow flood simulation 2.00 August 2015 modified for tsunami modelling in SETAS Author : E Rigby (0437 250 500) ted.rigby@rienco.com.au -------------------------------------------------------------------------------- """ ################################################################################# # # # IMPORT THE REQUIRED PYTHON MODULES # # # ################################################################################# # Standard python modules import os import sys from os import sep from os.path import join # Related ANUGA modules import anuga # Other python libraries import shapefile # the pyshp library - need to have shapefile.py accessible # Simulation specific imports import model_ini # Provides the simulation specific variables/files print model_ini.prefix+' >>>> Imported required modules for the model_data assembly process ' ################################################################################# # # # ASSEMBLE THE BACKDROP COASTAL POLYGON # # # ################################################################################# print model_ini.prefix+' >>>> Commencing data assembly for this simulation ' # load the coasline poly as backdrop for check plots if check true if model_ini.model_checks : sf = shapefile.Reader(model_ini.coastal_poly_pathname) # power up pyshp reader and get shape geom shapes =sf.shapes() print model_ini.prefix+' >>>> Reading in the backdrop coastline for check plots if checks true' print model_ini.prefix+' >>>> %i polygons in coastal shapefile ' %(len(shapes)) coastal_polygons =[] for i in range(len(shapes)) : # loop for each poly in shp file coastal_polygon_points = [] coastal_polygon_points = shapes[i].points # list of of xy coords for poly less closing point coastal_polygons.append(coastal_polygon_points) # append poly points list to polygons list if model_ini.model_verbose : for i in range(len(shapes)) : print model_ini.prefix+' ----- dump of first 2 coords of each coastal poly (i,[x,y]) - ', i,coastal_polygons[i][0:2] poly_window = anuga.plot_polygons(coastal_polygons, style = None, figname=join(model_ini.checks_dir,'Tasmanian Coastline'), label='Tasmanian Coastline') ################################################################################ # # # ASSEMBLE THE BOUNDING POLYGON & TAGGED BOUNDARY STRINGS # # # ################################################################################ #------------------------------------------------------------------------------ # Define bounding polygon and default res to be used later in model_build.py to # build the mesh and boundary strings with their tags # # OUTPUTS OF THIS CODE BLOCK ARE; # bounding_polygon[] -- The bounding polygon within which the mesh # will be created (from the bounding_poly shapefile) # boundary_strings{} -- The tagged boundary strings dictionary # #------------------------------------------------------------------------------ ##### Read in bounding poly and check plot if checks set to true ##### Note that the GA sts BC data contains boundary points in a counterclockwise direction ##### around the bounding poly starting from the most southerly ocean point. ##### As a consequence any attempt to create the bounding poly as a conventional GIS (shp) ##### file fails as the shape file automatically reorders in a clockwise manner around the bounding poly ##### as that is the required order for an outer ring/boundary in GIS. (counter clockwise denotes ##### an inner ring/hole poly). Must therefor read the bounding poly from a csv points list directly ##### into Anuga so that the non GIS standard of the sts oredering can be maintained. print model_ini.prefix+' >>>> Reading in the bounding polygon from ', model_ini.bounding_poly_pathname # single poly # this deleted untill/unless problem with sts order corrected. #sf = shapefile.Reader(model_ini.bounding_poly_pathname) # power up pyshp reader #shapes =sf.shapes() # read basic geometry of shapefile #if len(shapes) > 1 : # check only single poly read in # print model_ini.prefix+' ERROR : can only be one poly for bounding poly, i% read in ' %(len(shapes)) # sys.exit(0) #bounding_polygon = shapes[0].points[0:-1] # list of tuples of xy coords less closing point # bounding poly read direct from csv created from sum sts bdry points and land bdry points. # this creates an Anuga poly in the same order as the sts starting from the south ocean bdry # running north and then back down the land boundary in a non GS standad counterclockwise loop. bounding_polygon = anuga.read_polygon(model_ini.bounding_poly_pathname) if model_ini.model_verbose : print model_ini.prefix+' ----- There were %i unique points read in for the bounding polygon ' %(len(bounding_polygon)) print model_ini.prefix+' ----- dump of coords of the bounding poly [x,y]) - ', bounding_polygon if model_ini.model_checks : print model_ini.prefix+ ' >>>> Check plotting bounding polygon' # plot it out in data directory plot_polys = [] plot_polys = coastal_polygons plot_polys.append(bounding_polygon) poly_window = anuga.plot_polygons(plot_polys, figname=join(model_ini.checks_dir,'bounding_polygon'), label='Bounding Polygon') print model_ini.prefix+ ' >>>> Area of bounding polygon in Km2', anuga.polygon_area(bounding_polygon)/1000000.0 def_triangle_side = ((model_ini.default_res)* 2.0)**0.5 print model_ini.prefix+ ' >>>> Default resolution within bounding polygon is ',model_ini.default_res,' m2 ie ',def_triangle_side,' m triangle side' ##### Assign the required tagged boundary segments (verticies are those of bounding_poly) print model_ini.prefix+ ' >>>> Assembling the Anuga tagged boundary strings' # Note ANUGA will auto assign 'exterior' to any remaining segments of the bounding_poly # vertex numbers/segnents are derived from display of the bounding poly vertex numbers # Three explict and one implicit ('exterior') boundary strings are used to reflect the different BCs in the model # ANUGA will auto assign 'exterior' to any segments of the bounding_poly not assigned here. # BCs will be assigned to these boundary types in model_build.py # Note there are 56 gauges in the SETAS sts boundary file so 55 ocean_segments # the sts is indexed from the south to the north creating a non (GIS) # standard (counter clockwise) numbering sequence around the Anuga bounding_polygon # # top # viz 56*----<----*55 # | | # V ^ # Land | | Ocean (54 segments) # 64*----->---*0 # bottom top_segments = [] bottom_segments = [] ocean_segments = [] top_segments = [55] bottom_segments = [64] ocean_segments = range(55) # 0,1,,,54 boundary_strings = {'top' : top_segments, 'bottom' : bottom_segments, 'ocean' : ocean_segments} #if model_ini.model_checks : # need to sort out how to best check plot dictionary segments from vertices????ER?? # for key in boundary_strings : # plot_filename = model_ini.checks_dir +sep +'Model_Taggedbdry_Strings'+str(key) # bdry_segs_lst = boundary_strings[key] # plot_poly=[] # for i in range(len(bdry_segs_lst)) : # start_vertex = int(bdry_segs_lst[i]) # if start_vertex == len(bounding_polygon)-1 : # end_vertex = 0 # else : # end_vertex = start_vertex+1 # plot_poly.append(bounding_polygon[start_vertex]) # plot_poly.append(bounding_polygon[end_vertex]) # if model_ini.model_verbose : # print model_ini.prefix+ '----- plot bounding poly BC segments ' , key , plot_poly # poly_window = anuga.plot_polygons(plot_poly, style=None,figname=plot_filename) print model_ini.prefix+ ' >>>> %i tagged boundary strings assigned by user to model' %(len(boundary_strings)) if model_ini.model_verbose : for key in boundary_strings : print model_ini.prefix+ ' ----- tagged boundary ', key, boundary_strings[key] ################################################################################# # # # ASSEMBLE INTERIOR RESOLUTION POLYGONS # # # ################################################################################# #------------------------------------------------------------------------------- # Define Interior region polygons for later use in model_run.py to define zones # of different mesh resolution # # OUTPUTS OF THIS CODE BLOCK ARE; # interior_res_regions[]-- List of interior res regions containing [res_poly, res] # to be used in model_run to create the mesh #------------------------------------------------------------------------------- # Note coarse_polygon not needed as set as default for bounding_polgon print model_ini.prefix+ ' >>>> Assembling any inner region polygons and setting associated resolutions' # Interior_res_regions and assoc resolutions list are created here for use in model_build.py # order so coarsest interior res is first and finest last as each overstamps when assigned # note res polys must be simple, single part with vertices spaced out further than res. # A later res polys can be used when it lies fully inside a prior loaded poly but poly # boundaries canot touch or intersect. There should also be sufficient space between the # inner and outer polys for decently shaped triangles to form. # # the foillowing copde can probably be compacted into some sort of loop so can have variable # numbers of int polys ???????????ER interior_res_polygons = [] interior_res_regions = [] print model_ini.prefix+' >>>>> ',"about to load internal res1 medium ocean10to30m polygons ", #model_ini.int1_res_poly_pathname int1_polygons = [] sf = shapefile.Reader(model_ini.int1_res_poly_pathname) # power up pyshp reader shapes =sf.shapes() region_data = [] # read basic geometry of shapefile for i in range(len(shapes)) : # loop for each poly in shp file int1_polygon_points = [] int1_polygon_points = shapes[i].points # list of polys xy coords less closing point int1_polygons.append(int1_polygon_points) # append the points for that poly interior_res_polygons.append(int1_polygon_points) # keep polygons seperate for check plots region_data=[int1_polygon_points,model_ini.int1_res] # Add the int res data for this poly interior_res_regions.append(region_data) # Append these region_data polys to int_res_regions if model_ini.model_verbose : print model_ini.prefix+' ----- ', model_ini.int1_res_poly_pathname,' read in containing ', len(shapes), ' polygons' print model_ini.prefix+' >>>>> ',"about to load internal res2 fine coastal strip polygons ", model_ini.int2_res_poly_pathname int2_polygons = [] sf = shapefile.Reader(model_ini.int2_res_poly_pathname) shapes =sf.shapes() region_data = [] for i in range(len(shapes)) : int2_polygon_points = [] int2_polygon_points = shapes[i].points int2_polygons.append(int2_polygon_points) interior_res_polygons.append(int2_polygon_points) region_data = [int2_polygon_points, model_ini.int2_res] interior_res_regions.append(region_data) if model_ini.model_verbose : print model_ini.prefix+' ----- ', model_ini.int2_res_poly_pathname,' read in containing ', len(shapes), ' polygons' print model_ini.prefix+' >>>>> ',"about to load internal res3 coarse coastal land above 10m polygons",model_ini.int3_res_poly_pathname int3_polygons = [] sf = shapefile.Reader(model_ini.int3_res_poly_pathname) shapes =sf.shapes() region_data = [] for i in range(len(shapes)) : int3_polygon_points = [] int3_polygon_points = shapes[i].points int3_polygons.append(int3_polygon_points) interior_res_polygons.append(int3_polygon_points) region_data = [int3_polygon_points, model_ini.int3_res] interior_res_regions.append(region_data) if model_ini.model_verbose : print model_ini.prefix+' ----- ', model_ini.int3_res_poly_pathname,' read in containing ', len(shapes), ' polygons' print model_ini.prefix+' >>>>> ',"about to load internal res4 vfine plot area polygons ", model_ini.int4_res_poly_pathname int4_polygons = [] sf = shapefile.Reader(model_ini.int4_res_poly_pathname) shapes =sf.shapes() region_data = [] for i in range(len(shapes)) : int4_polygon_points = [] int4_polygon_points = shapes[i].points int4_polygons.append(int4_polygon_points) interior_res_polygons.append(int4_polygon_points) region_data = [int4_polygon_points, model_ini.int4_res] interior_res_regions.append(region_data) if model_ini.model_verbose : print model_ini.prefix+' ----- ', model_ini.int4_res_poly_pathname,' read in containing ', len(shapes), ' polygons' print model_ini.prefix+' >>>>> ',"about to load internal res5 xfine study area polygons ", model_ini.int5_res_poly_pathname int5_polygons = [] sf = shapefile.Reader(model_ini.int5_res_poly_pathname) shapes =sf.shapes() region_data = [] for i in range(len(shapes)) : int5_polygon_points = [] int5_polygon_points = shapes[i].points int5_polygons.append(int5_polygon_points) interior_res_polygons.append(int5_polygon_points) region_data = [int5_polygon_points, model_ini.int5_res] interior_res_regions.append(region_data) if model_ini.model_verbose : print model_ini.prefix+' ----- ', model_ini.int5_res_poly_pathname,' read in containing ', len(shapes), ' polygons' # Do some final sanity checks on the poly data read in if (len(interior_res_regions) != len(interior_res_polygons)) : print model_ini.prefix+ ' ERROR - Length of int_respolys and int_res_regions not the same ', len(interior_res_polygons), len(interior_res_regions) sys.exit(0) print model_ini.prefix+ ' >>>> %i interior resolution polygons created and processed' %(len(interior_res_regions)) if model_ini.model_verbose : for i in range(len(interior_res_regions)) : print model_ini.prefix+' ----- dump of internal resolution m2 for poly ', i, interior_res_regions[i][1] for i in range(len(interior_res_regions)) : print model_ini.prefix+' ----- dump of coords for poly ', i, interior_res_polygons[i][0],'...',interior_res_polygons[i][-1] if model_ini.model_checks : print model_ini.prefix+ ' >>>> Check plotting interior resolution polygons' # plot it out in data directory plot_polys = [] plot_polys= interior_res_polygons # plot_polys.append(bounding_polygon) # seems to create a problem if appended???ER?? poly_window = anuga.plot_polygons(plot_polys, figname=join(model_ini.checks_dir,'Interior Resolution Polygons'), label='Interior Resolution Polygons') ################################################################################# # # # ASSEMBLE THE ELEVATION DATA # # # ################################################################################# #------------------------------------------------------------------------------ # Assign reference to file with elevation data for later use in model_build.py # the following assumes this elevation data file is stored with the other # datafiles in the model/scenario/event/data directory # # OUTPUTS OF THIS CODE BLOCK ARE; # elevation_pathname -- full path name of file containing ALS data to be used to set mesh elevations #------------------------------------------------------------------------------ # Note must be either 1 or 2 #elev_option = 1 # assign elevations to the mesh using the fast netcdf elevation dataset elev_option = 2 # assign elevations to the mesh using the original (slower) yxz/txt dataset if elev_option != 1 and elev_option != 2 : print model_ini.prefix+ ' >>>> elev_option must be 1 or 2 - correct and rerun - was ', elev_option # NOTES ##### # OPTION 1 (Preferred but limited to smaller ALS datasets) # This option converts the raw xyz ALS data into NETCDF (pts) format # The pts file is then used to assign mesh elevations in model_build.py. # While interpolating and assigning elevations to the mesh from the pts file, # this option loads all data into memory (no blocking) so is limited in # regard to the number of points it can assign to the mesh. # If there are more than about 6 million points in the ALS dataset # this option will fail with a memory over-run error. # In such a case the elevation data should be left in csv (xyz/txt) format. # # OPTION 2 (Leave elevation data in original (csv)(xyz) format and block read when elevating mesh) # if elev_option == 1 : # Convert xyz/txt data into NETCDF pts format to speed later use in setting mesh elevations print model_ini.prefix+' >>>> Elevation x,y,z scatterpoint dataset %s being used to create NETCDF pts file' %(model_ini.elevation_points_pathname) G = anuga.Geospatial_data(file_name = model_ini.elevation_points_pathname, max_read_lines = 500, load_file_now = True, verbose = model_ini.model_verbose) G.export_points_file(model_ini.elevation_points_pathname[:-4] + '.pts') # export geospatial data to NETCDF pts formatted file elevation_pathname = model_ini.elevation_points_pathname[:-4] + '.pts' # switch elevation_filename to now point to the NETCDF pts file print model_ini.prefix+ ' >>>> Mesh elevation data will be assigned from %s' %(elevation_pathname) else : # OPTION 2 (Leave elevation data in original (csv)(xyz) format and block read when elevating mesh in model_build.py) elevation_pathname = model_ini.elevation_points_pathname print model_ini.prefix+ ' >>>> Mesh elevation data will be assigned from %s' %(elevation_pathname) ################################################################################# # # # ASSEMBLE SURFACE ROUGHNESS POLYGONS AND ASSOCIATED DATA # # # ################################################################################# # This block would typically read in the roughness data from a shapefile but for # more complex roughness dtaasets the roughness will be directly assigned to the # mesh in model_build from an asc grid file print model_ini.prefix+ ' >>>> The surface roughness data will be assembled from ' , model_ini.mannings_grid_pathname ################################################################################ # # # ASSEMBLE THE SPATIO/TEMPORAL BOUNDARY CONDITIONS # # # ################################################################################ #------------------------------------------------------------------------------- # # Note : For this model the temporal conditions (stage) at the sea boundary # were preprocessed by GA in 2009 using build_boundaries.py (no longer available) # into a series of sts (binary netcdf) files for each tsunami event). # While this permits this newer model to be run with the same BC including various # fixed tidal variation via the optional mean_stage parameter, it may present # a problem when trying to combine a temporally varying tide to the tsunami. # ???????????????Needs further review????????????????????ER??? # #-------------------------------------------------------------------------------