"""Build the elevation data to run a storm surge inundation scenario for Mandurah, Western Australia, Australia. Input: elevation data from project.py Output: pts file stored in project.topographies_dir The run_model.py is reliant on the output of this script. """ #------------------------------------------------------------------------------ # Import necessary modules #------------------------------------------------------------------------------ # Standard modules import os from os.path import join # Related major packages from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf from anuga.shallow_water.data_manager import dem2pts from anuga.shallow_water.data_manager import start_screen_catcher from anuga.shallow_water.data_manager import copy_code_files from anuga.geospatial_data.geospatial_data import Geospatial_data # Application specific imports from setup_model import project # Definition of file names and polygons #------------------------------------------------------------------------------ # Copy scripts to time stamped output directory and capture screen # output to file #------------------------------------------------------------------------------ copy_code_files(project.output_build,__file__, os.path.dirname(project.__file__)+os.sep+\ project.__name__+'.py' ) start_screen_catcher(project.output_build) #------------------------------------------------------------------------------ # 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.bounding_polygon', project.bounding_polygon print 'project.combined_elevation_basename', project.combined_elevation_basename # Create Geospatial data from ASCII files geospatial_data = {} if not project.ascii_grid_filenames == []: for filename in project.ascii_grid_filenames: absolute_filename = join(project.topographies_folder, filename) convert_dem_from_ascii2netcdf(absolute_filename, basename_out=absolute_filename, use_cache=True, verbose=True) dem2pts(absolute_filename, use_cache=True, verbose=True) G_grid = Geospatial_data(file_name=absolute_filename+'.pts', verbose=True) print 'Clip geospatial object' geospatial_data[filename] = G_grid.clip(project.bounding_polygon) # Create Geospatial data from TXT files if not project.point_filenames == []: for filename in project.point_filenames: absolute_filename = join(project.topographies_folder, filename) G_points = Geospatial_data(file_name=absolute_filename, verbose=True) print 'Clip geospatial object' geospatial_data[filename] = G_points.clip(project.bounding_polygon) print 'Geospatial data objects: ', geospatial_data.keys() #------------------------------------------------------------------------------- # Combine, clip and export dataset #------------------------------------------------------------------------------- print 'Add geospatial objects' # except', project.offshore_name5 print 'project.elevation_clip_box', project.elevation_clip_box keylist = project.ascii_grid_filenames + project.point_filenames G1 = geospatial_data[keylist[0]].clip_outside(project.elevation_clip_box) G2 = geospatial_data[keylist[1]].clip(project.elevation_clip_box) G3 = geospatial_data[keylist[2]] G4 = geospatial_data[keylist[3]] G5 = geospatial_data[keylist[4]] G6 = geospatial_data[keylist[5]] G7 = geospatial_data[keylist[6]] G8 = geospatial_data[keylist[7]] G9 = geospatial_data[keylist[8]] G10 = geospatial_data[keylist[9]] #G11 = geospatial_data[keylist[10]] # G = None # for i, key in enumerate(keylist): # if key==keylist[0]: # G[i] = geospatial_data[key].clip_outside(project.elevation_clip_box) # elif key==keylist[1]: # G[i] = geospatial_data[key].clip(project.elevation_clip_box) # else # G+ = geospatial_data[key] #~ G_list = [G1, G2, G3, G4, G5, G6, G7, G8, G9, G10] #~ print 'G1', G1.attributes #~ print 'G2', G2.attributes #~ print 'G3', G3.attributes #~ print 'G4', G4.attributes #~ print 'G5', G5.attributes #~ print 'G6', G6.attributes #~ print 'G7', G7.attributes #~ print 'G8', G8.attributes #~ print 'G9', G9.attributes #~ print 'G10', G10.attributes G = G1 + G2 + G3 + G4 + G5 + G6 + G7 + G8 + G9 + G10 #+ G11 # G = None # for key in geospatial_data: # if key == project.ascii_grid_filenames[0]: # G == geospatial_data[key].clip_outside(project.elevation_clip_box) # elif key == project.ascii_grid_filenames[1]: # G = geospatial_data[key].clip(project.elevation_clip_box) # continue # G += geospatial_data[key] print 'Export combined DEM file' G.export_points_file(project.combined_elevation + '.pts') print 'Do txt version too' # Use for comparision in ARC #~ for g in G_list: #~ for i in keylist: #~ g.export_points_file(join(project.topographies_folder, (str(keylist[i]) +'_export.txt'))) G.export_points_file(project.combined_elevation + '.txt') print 'Export individual text files' #for i in xrange(len(keylist)): # G[i+1].export_points_file(join(project.topographies_folder, keylist[i] +'_export.txt')) G1.export_points_file(join(project.topographies_folder, keylist[0] +'_export.txt')) G2.export_points_file(join(project.topographies_folder, keylist[1] +'_export.txt')) G3.export_points_file(join(project.topographies_folder, keylist[2] +'_export.txt')) G4.export_points_file(join(project.topographies_folder, keylist[3] +'_export.txt')) G5.export_points_file(join(project.topographies_folder, keylist[4] +'_export.txt')) G6.export_points_file(join(project.topographies_folder, keylist[5] +'_export.txt')) G7.export_points_file(join(project.topographies_folder, keylist[6] +'_export.txt')) G8.export_points_file(join(project.topographies_folder, keylist[7] +'_export.txt')) G9.export_points_file(join(project.topographies_folder, keylist[8] +'_export.txt')) G10.export_points_file(join(project.topographies_folder, keylist[9] +'_export.txt')) #G11.export_points_file(join(project.topographies_folder, keylist[10] +'_export.txt'))