[7613] | 1 | """Build the elevation data to run a storm surge inundation scenario |
---|
| 2 | for Mandurah, Western Australia, Australia. |
---|
| 3 | |
---|
| 4 | Input: elevation data from project.py |
---|
| 5 | Output: pts file stored in project.topographies_dir |
---|
| 6 | The run_model.py is reliant on the output of this script. |
---|
| 7 | |
---|
| 8 | """ |
---|
| 9 | |
---|
| 10 | #------------------------------------------------------------------------------ |
---|
| 11 | # Import necessary modules |
---|
| 12 | #------------------------------------------------------------------------------ |
---|
| 13 | |
---|
| 14 | # Standard modules |
---|
| 15 | import os |
---|
| 16 | from os.path import join |
---|
| 17 | |
---|
| 18 | # Related major packages |
---|
| 19 | from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf |
---|
| 20 | from anuga.shallow_water.data_manager import dem2pts |
---|
| 21 | from anuga.shallow_water.data_manager import start_screen_catcher |
---|
| 22 | from anuga.shallow_water.data_manager import copy_code_files |
---|
| 23 | from anuga.geospatial_data.geospatial_data import Geospatial_data |
---|
| 24 | |
---|
| 25 | # Application specific imports |
---|
| 26 | from setup_model import project # Definition of file names and polygons |
---|
| 27 | |
---|
| 28 | |
---|
| 29 | #------------------------------------------------------------------------------ |
---|
| 30 | # Copy scripts to time stamped output directory and capture screen |
---|
| 31 | # output to file |
---|
| 32 | #------------------------------------------------------------------------------ |
---|
| 33 | copy_code_files(project.output_build,__file__, |
---|
| 34 | os.path.dirname(project.__file__)+os.sep+\ |
---|
| 35 | project.__name__+'.py' ) |
---|
| 36 | start_screen_catcher(project.output_build) |
---|
| 37 | |
---|
| 38 | #------------------------------------------------------------------------------ |
---|
| 39 | # Preparation of topographic data |
---|
| 40 | # |
---|
| 41 | # Convert ASC 2 DEM 2 PTS using source data and store result in source data |
---|
| 42 | # Do for coarse and fine data |
---|
| 43 | # Fine pts file to be clipped to area of interest |
---|
| 44 | #------------------------------------------------------------------------------ |
---|
| 45 | |
---|
| 46 | print 'project.bounding_polygon', project.bounding_polygon |
---|
| 47 | print 'project.combined_elevation_basename', project.combined_elevation_basename |
---|
| 48 | |
---|
| 49 | # Create Geospatial data from ASCII files |
---|
| 50 | geospatial_data = {} |
---|
| 51 | if not project.ascii_grid_filenames == []: |
---|
| 52 | for filename in project.ascii_grid_filenames: |
---|
| 53 | absolute_filename = join(project.topographies_folder, filename) |
---|
| 54 | convert_dem_from_ascii2netcdf(absolute_filename, |
---|
| 55 | basename_out=absolute_filename, |
---|
| 56 | use_cache=True, |
---|
| 57 | verbose=True) |
---|
| 58 | dem2pts(absolute_filename, use_cache=True, verbose=True) |
---|
| 59 | |
---|
| 60 | G_grid = Geospatial_data(file_name=absolute_filename+'.pts', |
---|
| 61 | verbose=True) |
---|
| 62 | print 'Clip geospatial object' |
---|
| 63 | geospatial_data[filename] = G_grid.clip(project.bounding_polygon) |
---|
| 64 | |
---|
| 65 | # Create Geospatial data from TXT files |
---|
| 66 | if not project.point_filenames == []: |
---|
| 67 | for filename in project.point_filenames: |
---|
| 68 | absolute_filename = join(project.topographies_folder, filename) |
---|
| 69 | G_points = Geospatial_data(file_name=absolute_filename, |
---|
| 70 | verbose=True) |
---|
| 71 | print 'Clip geospatial object' |
---|
| 72 | geospatial_data[filename] = G_points.clip(project.bounding_polygon) |
---|
| 73 | |
---|
| 74 | print 'Geospatial data objects: ', geospatial_data.keys() |
---|
| 75 | |
---|
| 76 | #------------------------------------------------------------------------------- |
---|
| 77 | # Combine, clip and export dataset |
---|
| 78 | #------------------------------------------------------------------------------- |
---|
| 79 | |
---|
| 80 | print 'Add geospatial objects' # except', project.offshore_name5 |
---|
| 81 | |
---|
| 82 | print 'project.elevation_clip_box', project.elevation_clip_box |
---|
| 83 | |
---|
| 84 | keylist = project.ascii_grid_filenames + project.point_filenames |
---|
| 85 | |
---|
[7647] | 86 | G1 = geospatial_data['swwa_10m_IC'].clip_outside(project.elevation_clip_box) |
---|
| 87 | G2 = geospatial_data['bunbury_5m_IC'].clip(project.elevation_clip_box) |
---|
[7613] | 88 | #~ G2 = geospatial_data['bunbury_nth_a'].clip(project.elevation_clip_box) |
---|
| 89 | #~ G3 = geospatial_data['bunbury_nth_b'].clip(project.elevation_clip_box) |
---|
| 90 | #~ G4 = geospatial_data['bunbury_sth_a'].clip(project.elevation_clip_box) |
---|
| 91 | #~ G5 = geospatial_data['bunbury_sth_b'].clip(project.elevation_clip_box) |
---|
| 92 | G3 = geospatial_data['DPI.txt'] |
---|
| 93 | G4 = geospatial_data['Busselton_Chart_Clip_ss.txt'] |
---|
| 94 | G5 = geospatial_data['Busselton_NavyFinal_Clip_ss.txt'] |
---|
| 95 | G6 = geospatial_data['Leschenault_TIN.txt'] |
---|
[7647] | 96 | G7 = geospatial_data['Leschenault_Inlet_TIN.txt'] |
---|
| 97 | G8 = geospatial_data['DPI5U1A02_01a_edited.txt'] |
---|
| 98 | G9 = geospatial_data['DPI5U1A02_01b_edited.txt'] |
---|
| 99 | G10 = geospatial_data['DPI5U1A02_01c_edited.txt'] |
---|
| 100 | G11 = geospatial_data['DPI5U1A02_01d_edited.txt'] |
---|
| 101 | G12 = geospatial_data['DPI5U1A02_01e_edited.txt'] |
---|
[7613] | 102 | |
---|
| 103 | #~ G_list = [G1, G2, G3, G4, G5, G6, G7, G8, G9, G10] |
---|
| 104 | |
---|
| 105 | #~ print 'G1', G1.attributes |
---|
| 106 | #~ print 'G2', G2.attributes |
---|
| 107 | #~ print 'G3', G3.attributes |
---|
| 108 | #~ print 'G4', G4.attributes |
---|
| 109 | #~ print 'G5', G5.attributes |
---|
| 110 | #~ print 'G6', G6.attributes |
---|
| 111 | #~ print 'G7', G7.attributes |
---|
| 112 | #~ print 'G8', G8.attributes |
---|
| 113 | #~ print 'G9', G9.attributes |
---|
| 114 | #~ print 'G10', G10.attributes |
---|
| 115 | |
---|
[7647] | 116 | G = G1 + G2 + G3 + G4 + G5 + G6 + G7 + G8 + G9 + G10 + G11 + G12 |
---|
[7613] | 117 | |
---|
| 118 | # G = None |
---|
| 119 | # for key in geospatial_data: |
---|
| 120 | # if key == project.ascii_grid_filenames[0]: |
---|
| 121 | # G == geospatial_data[key].clip_outside(project.elevation_clip_box) |
---|
| 122 | # elif key == project.ascii_grid_filenames[1]: |
---|
| 123 | # G = geospatial_data[key].clip(project.elevation_clip_box) |
---|
| 124 | # continue |
---|
| 125 | |
---|
| 126 | # G += geospatial_data[key] |
---|
| 127 | |
---|
| 128 | print 'Export combined DEM file' |
---|
| 129 | G.export_points_file(project.combined_elevation + '.pts') |
---|
| 130 | |
---|
| 131 | print 'Do txt version too' |
---|
| 132 | # Use for comparision in ARC |
---|
| 133 | #~ for g in G_list: |
---|
| 134 | #~ for i in keylist: |
---|
| 135 | #~ g.export_points_file(join(project.topographies_folder, (str(keylist[i]) +'_export.txt'))) |
---|
| 136 | G.export_points_file(project.combined_elevation + '.txt') |
---|
| 137 | G1.export_points_file(join(project.topographies_folder, keylist[0] +'_export.txt')) |
---|
| 138 | G2.export_points_file(join(project.topographies_folder, keylist[1] +'_export.txt')) |
---|
| 139 | G3.export_points_file(join(project.topographies_folder, keylist[2] +'_export.txt')) |
---|
| 140 | G4.export_points_file(join(project.topographies_folder, keylist[3] +'_export.txt')) |
---|
| 141 | G5.export_points_file(join(project.topographies_folder, keylist[4] +'_export.txt')) |
---|
| 142 | G6.export_points_file(join(project.topographies_folder, keylist[5] +'_export.txt')) |
---|
| 143 | G7.export_points_file(join(project.topographies_folder, keylist[6] +'_export.txt')) |
---|
| 144 | G8.export_points_file(join(project.topographies_folder, keylist[7] +'_export.txt')) |
---|
| 145 | G9.export_points_file(join(project.topographies_folder, keylist[8] +'_export.txt')) |
---|
| 146 | G10.export_points_file(join(project.topographies_folder, keylist[9] +'_export.txt')) |
---|
[7647] | 147 | G11.export_points_file(join(project.topographies_folder, keylist[10] +'_export.txt')) |
---|
| 148 | G12.export_points_file(join(project.topographies_folder, keylist[11] +'_export.txt')) |
---|