"""Script for running tsunami inundation scenario for Perth, WA, Australia. Source data such as elevation and boundary data is assumed to be available in directories specified by project_clip.py The output sww file is stored in project_clip.output_time_dir The scenario is defined by a triangular mesh created from project_clip.polygon, the elevation data and a simulated submarine landslide. Ole Nielsen and Duncan Gray, GA - 2005 and Jane Sexton, Nick Bartzis, GA - 2006 """ #------------------------------------------------------------------------------ # Import necessary modules #------------------------------------------------------------------------------ # Standard modules from os import sep from os.path import dirname, basename from os import mkdir, access, F_OK from shutil import copy import time import sys # Related major packages from anuga.shallow_water import Domain #from anuga.shallow_water import Dirichlet_boundary #from anuga.shallow_water import File_boundary #from anuga.shallow_water import Reflective_boundary from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf, dem2pts #from anuga.pmesh.mesh_interface import create_mesh_from_regions from anuga.geospatial_data.geospatial_data import * from anuga.shallow_water.data_manager import start_screen_catcher, copy_code_files # Application specific imports import project_clip # Definition of file names and polygons #------------------------------------------------------------------------------ # Copy scripts to time stamped output directory and capture screen # output to file #------------------------------------------------------------------------------ copy_code_files(project_clip.output_build_time_dir,__file__, dirname(project_clip.__file__)+sep+ project_clip.__name__+'.py' ) start_screen_catcher(project_clip.output_build_time_dir) print 'USER: ', project_clip.user #------------------------------------------------------------------------------- # Preparation of topographic data# Convert ASC 2 DEM 2 PTS using source data and store result in source data # Do for coarse and fine data # Fine pts file to be clipped to area of interest #------------------------------------------------------------------------------- print"project_clip.combined_dir_name",project_clip.combined_dir_name # topography directory filenames total=project_clip.combined_dir_name print'create Geospatial data1 objects from topographies',total + '.txt' G1 = Geospatial_data(file_name = total + '.txt') print'add all geospatial objects' Gt= G1.clip_outside(Geospatial_data(project_clip.poly_west_reef)).clip_outside(Geospatial_data(project_clip.poly_west)).clip_outside(Geospatial_data(project_clip.poly_centerKolonia)).clip_outside(Geospatial_data(project_clip.poly_east)) #print'clip combined geospatial object by bounding polygon' #G_clipped = G.clip(project_clip_urs.poly_all) #FIXME: add a clip function to pts #print'shape of clipped data', G_clipped.get_data_points().shape print'export combined DEM file' if access(project_clip.topographies_dir,F_OK) == 0: mkdir (project_clip.topographies_dir) print'export',project_clip.combined_dir_name+ '.txt' Gt.export_points_file(project_clip.combined_dir_name1+ '.txt') ''' print'project_clip.combined_dir_name + 1.xya',project_clip.combined_dir_name + '1.xya' G_all=Geospatial_data(file_name = project_clip.combined_dir_name + '1.xya') print'split' G_all_1, G_all_2 = G_all.split(.10) print'export 1' G_all_1.export_points_file(project_clip.combined_dir_name+'_small1' + '.xya') print'export 2' G_all_2.export_points_file(project_clip.combined_dir_name+'_other1' + '.xya') #------------------------------------------------------------------------- # Convert URS to SWW file for boundary conditions #------------------------------------------------------------------------- print 'starting to create boundary conditions' boundaries_in_dir_name = project_clip.boundaries_in_dir_name from anuga.shallow_water.data_manager import urs2sww print 'minlat=project_clip.north_boundary, maxlat=project_clip.south_boundary',project_clip.north_boundary, project_clip.south_boundary print 'minlon= project_clip.west_boundary, maxlon=project_clip.east_boundary',project_clip.west_boundary, project_clip.east_boundary #import sys; sys.exit() #if access(project_clip.boundaries_dir,F_OK) == 0: # mkdir (project_clip.boundaries_dir) from caching import cache cache(urs2sww, (boundaries_in_dir_name, project_clip.boundaries_dir_name1), {'verbose': True, 'minlat': project_clip.south_boundary, 'maxlat': project_clip.north_boundary, 'minlon': project_clip.west_boundary, 'maxlon': project_clip.east_boundary, # 'minlat': project_clip.south, # 'maxlat': project_clip.north, # 'minlon': project_clip.west, # 'maxlon': project_clip.east, 'mint': 0, 'maxt': 40000, # 'origin': domain.geo_reference.get_origin(), 'mean_stage': project_clip.tide, # 'zscale': 1, #Enhance tsunami 'fail_on_NaN': False}, verbose = True, ) # dependencies = source_dir + project_clip.boundary_basename + '.sww') '''