"""Script for running a tsunami inundation scenario for Cairns, QLD Australia. Source data such as elevation and boundary data is assumed to be available in directories specified by project.py The output sww file is stored in directory named after the scenario, i.e slide or fixed_wave. The scenario is defined by a triangular mesh created from project.polygon, the elevation data and a tsunami wave generated by a submarine mass failure. Ole Nielsen and Duncan Gray, GA - 2005 and Jane Sexton and Nick Bartzis, GA - 2006 """ def convert_arcgis_latlon_list_to_utm(points): #Used because arc gis produced csv files put lat lon in #reverse order to those accpeted by convert_latlon_to_utm() from anuga.coordinate_transforms.geo_reference import Geo_reference from anuga.coordinate_transforms.redfearn import redfearn old_geo = Geo_reference() utm_points = [] for point in points: zone, easting, northing = redfearn(float(point[1]),float(point[0])) new_geo = Geo_reference(zone) old_geo.reconcile_zones(new_geo) utm_points.append([easting,northing]) return utm_points, old_geo.get_zone() #------------------------------------------------------------------------------ # Import necessary modules #------------------------------------------------------------------------------ # Standard modules import os import time import sys # Related major packages from anuga.shallow_water import Domain from anuga.shallow_water import Reflective_boundary from anuga.shallow_water import Dirichlet_boundary from anuga.shallow_water import Time_boundary from anuga.shallow_water import File_boundary from anuga.pmesh.mesh_interface import create_mesh_from_regions from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf, ferret2sww from anuga.shallow_water.data_manager import dem2pts from anuga.coordinate_transforms.redfearn import convert_from_latlon_to_utm #from anuga.fit_interpolate.fit import fit_to_mesh_file (does not exist) # Application specific imports import project # Definition of file names and polygons #------------------------------------------------------------------------------ # Define scenario as either slide or fixed_wave. #------------------------------------------------------------------------------ #scenario = 'coseismic' scenario = 'fixed_wave' if os.access(scenario, os.F_OK) == 0: os.mkdir(scenario) basename = scenario + 'source' #------------------------------------------------------------------------------ # Preparation of topographic data # Convert ASC 2 DEM 2 PTS using source data and store result in source data #------------------------------------------------------------------------------ # Filenames dem_name = 'boxingday' # Create DEM from asc data #convert_dem_from_ascii2netcdf(dem_name, use_cache=True, verbose=True) # Create pts file for onshore DEM #dem2pts(dem_name, use_cache=True, verbose=True) #------------------------------------------------------------------------------ # Create the triangular mesh based on overall clipping polygon with a tagged # boundary and interior regions defined in project.py along with # resolutions (maximal area of per triangle) for each polygon #------------------------------------------------------------------------------ remainder_res = 10000000 islands_res = 100000 cairns_res = 100000 shallow_res = 500000 #interior_regions = [[project.poly_cairns, cairns_res]#, # [project.poly_island0, islands_res], # [project.poly_island1, islands_res], # [project.poly_island2, islands_res], # [project.poly_island3, islands_res], # [project.poly_shallow, shallow_res]] # filenames meshname = project.meshname + '.tsh' mesh_elevname = project.mesh_elevname + '.tsh' bounding_polygon, zone = convert_arcgis_latlon_list_to_utm(project.bounding_polygon_latlon) from caching import cache print 'start create mesh from regions' if scenario == 'coseismic': _ = cache(create_mesh_from_regions, bounding_polygon, {'boundary_tags': {'one': [0], 'two': [1], 'three': [2], 'four': [3], 'five': [4], 'six': [5], 'seven': [6], 'eight': [7], 'nine': [8], 'ten': [9], 'eleven': [10], 'twelve': [11], 'thirteen': [12], 'fourteen': [13]}, 'maximum_triangle_area': 5000000, 'filename': meshname},#, #'interior_regions': interior_regions}, verbose = True #, evaluate=True ) if scenario == 'fixed_wave': south = 7.40911081272 north = 8.71484358635 east = 99.1683687224 west = 97.513856322 points = [[south,west],[south,east],[north,east],[north,west]] bounding_polygon,zone=convert_from_latlon_to_utm(points) _ = cache(create_mesh_from_regions, bounding_polygon, {'boundary_tags': {'ocean_west': [0], 'bottom': [1], 'onshore': [2], 'top': [3]}, 'maximum_triangle_area': 500000000, 'filename': meshname},#, #'interior_regions': interior_regions}, verbose = True #, evaluate=True ) ## cache(fit_to_mesh_file,(meshname, ## project.combined_dir_name + '.pts', ## mesh_elevname), ## #{'verbose': True}, ## verbose = False ## ) #------------------------------------------------------------------------------ # Setup computational domain #------------------------------------------------------------------------------ domain = Domain(meshname, use_cache=False, verbose=True) print 'Number of triangles = ', len(domain) print 'The extent is ', domain.get_extent() print domain.statistics() domain.set_name(basename) domain.set_datadir(scenario) domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum']) domain.set_minimum_storable_height(0.01) #------------------------------------------------------------------------------ # Setup initial conditions #------------------------------------------------------------------------------ def elevation(x,y): return -10.0 tide = 0.0 domain.set_quantity('stage', tide) domain.set_quantity('friction', 0.0) if scenario == 'fixed_wave': # test with coarser bathymetry domain.set_quantity('elevation', #filename = 'thaicoas_9.pts', elevation, use_cache=True, verbose=True, alpha=0.1) if scenario == 'coseismic': domain.set_quantity('elevation', filename=project.combined_dir_name + '.pts', use_cache=True, verbose=True, alpha=0.1) #------------------------------------------------------------------------------ # Setup boundary conditions #------------------------------------------------------------------------------ print 'start ferret2sww' south = project.south north = project.north west = project.west east = project.east """ tide = 0.0 # Note only need to do when an SWW file for the MOST boundary doesn't exist cache(ferret2sww, (project.boundary_most_in, project.boundary_most_out), {'verbose': True, 'minlat': south, 'maxlat': north, 'minlon': west, 'maxlon': east, # 'origin': project.mesh_origin, 'origin': domain.geo_reference.get_origin(), 'mean_stage': tide, 'zscale': 10, #Enhance tsunami (Does this affect result) 'fail_on_NaN': False, 'inverted_bathymetry': True}, #evaluate = True, verbose = True)#, #dependencies = 'most_results' + '.sww') """ print 'Available boundary tags', domain.get_boundary_tags() #Bf = File_boundary('most_results'+'.sww', # domain, verbose = True) Br = Reflective_boundary(domain) Bd = Dirichlet_boundary([tide,0,0]) Bw = Time_boundary(domain = domain, f=lambda t: [(60