[5810] | 1 | """ |
---|
| 2 | Script for building the elevation data to run tsunami inundation scenario |
---|
| 3 | for Broome, WA, Australia. |
---|
| 4 | |
---|
| 5 | Input: elevation data from project.py |
---|
| 6 | Output: pts file stored in project.topographies_dir |
---|
| 7 | The run_broome.py is reliant on the output of this script. |
---|
| 8 | |
---|
| 9 | Ole Nielsen and Duncan Gray, GA - 2005 and Jane Sexton, Nick Bartzis, GA - 2006 |
---|
| 10 | """ |
---|
| 11 | |
---|
| 12 | #------------------------------------------------------------------------------ |
---|
| 13 | # Import necessary modules |
---|
| 14 | #------------------------------------------------------------------------------ |
---|
| 15 | |
---|
| 16 | # Standard modules |
---|
| 17 | from os import sep |
---|
| 18 | from os.path import dirname, basename |
---|
| 19 | from os import mkdir, access, F_OK |
---|
| 20 | from shutil import copy |
---|
| 21 | import time |
---|
| 22 | import sys |
---|
| 23 | |
---|
| 24 | # Related major packages |
---|
| 25 | from anuga.shallow_water import Domain |
---|
| 26 | from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf, dem2pts |
---|
| 27 | from anuga.geospatial_data.geospatial_data import * |
---|
| 28 | from anuga.shallow_water.data_manager import start_screen_catcher, copy_code_files |
---|
| 29 | |
---|
| 30 | # Application specific imports |
---|
| 31 | import project # Definition of file names and polygons |
---|
| 32 | |
---|
| 33 | #------------------------------------------------------------------------------ |
---|
| 34 | # Copy scripts to time stamped output directory and capture screen |
---|
| 35 | # output to file |
---|
| 36 | #------------------------------------------------------------------------------ |
---|
| 37 | |
---|
| 38 | copy_code_files(project.output_build_time_dir,__file__, |
---|
| 39 | dirname(project.__file__)+sep+ project.__name__+'.py' ) |
---|
| 40 | |
---|
| 41 | start_screen_catcher(project.output_build_time_dir) |
---|
| 42 | |
---|
| 43 | print 'USER: ', project.user |
---|
| 44 | |
---|
| 45 | #------------------------------------------------------------------------------- |
---|
| 46 | # Preparation of topographic data |
---|
| 47 | # |
---|
| 48 | # Convert ASC 2 DEM 2 PTS using source data and store result in source data |
---|
| 49 | # Do for coarse and fine data |
---|
| 50 | # Fine pts file to be clipped to area of interest |
---|
| 51 | #------------------------------------------------------------------------------- |
---|
| 52 | print"project.poly_all",project.poly_all |
---|
| 53 | print"project.combined_dir_name",project.combined_dir_name |
---|
| 54 | |
---|
| 55 | # input elevation directory filenames |
---|
| 56 | onshore_in_dir_name = project.onshore_in_dir_name |
---|
| 57 | |
---|
| 58 | coast_in_dir_name = project.coast_in_dir_name |
---|
| 59 | |
---|
| 60 | offshore_in_dir_name = project.offshore_in_dir_name |
---|
| 61 | offshore_in_dir_name1 = project.offshore_in_dir_name1 |
---|
| 62 | offshore_in_dir_name2 = project.offshore_in_dir_name2 |
---|
| 63 | offshore_in_dir_name3 = project.offshore_in_dir_name3 |
---|
| 64 | offshore_in_dir_name4 = project.offshore_in_dir_name4 |
---|
| 65 | offshore_in_dir_name5 = project.offshore_in_dir_name5 |
---|
| 66 | offshore_in_dir_name6 = project.offshore_in_dir_name6 |
---|
| 67 | |
---|
| 68 | # output elevation directory filenames |
---|
| 69 | onshore_dir_name = project.onshore_dir_name |
---|
| 70 | |
---|
| 71 | coast_dir_name = project.coast_dir_name |
---|
| 72 | |
---|
| 73 | offshore_dir_name = project.offshore_dir_name |
---|
| 74 | offshore_dir_name1 = project.offshore_dir_name1 |
---|
| 75 | offshore_dir_name2 = project.offshore_dir_name2 |
---|
| 76 | offshore_dir_name3 = project.offshore_dir_name3 |
---|
| 77 | offshore_dir_name4 = project.offshore_dir_name4 |
---|
| 78 | offshore_dir_name5 = project.offshore_dir_name5 |
---|
| 79 | offshore_dir_name6 = project.offshore_dir_name6 |
---|
| 80 | |
---|
| 81 | # creates DEM from asc data |
---|
| 82 | print "creates DEMs from asc data" |
---|
| 83 | convert_dem_from_ascii2netcdf(onshore_in_dir_name, basename_out=onshore_dir_name, use_cache=True, verbose=True) |
---|
| 84 | |
---|
| 85 | convert_dem_from_ascii2netcdf(offshore_in_dir_name1, basename_out=offshore_dir_name1, use_cache=True, verbose=True) |
---|
| 86 | convert_dem_from_ascii2netcdf(offshore_in_dir_name2, basename_out=offshore_dir_name2, use_cache=True, verbose=True) |
---|
| 87 | convert_dem_from_ascii2netcdf(offshore_in_dir_name3, basename_out=offshore_dir_name3, use_cache=True, verbose=True) |
---|
| 88 | convert_dem_from_ascii2netcdf(offshore_in_dir_name4, basename_out=offshore_dir_name4, use_cache=True, verbose=True) |
---|
| 89 | convert_dem_from_ascii2netcdf(offshore_in_dir_name5, basename_out=offshore_dir_name5, use_cache=True, verbose=True) |
---|
| 90 | convert_dem_from_ascii2netcdf(offshore_in_dir_name6, basename_out=offshore_dir_name6, use_cache=True, verbose=True) |
---|
| 91 | |
---|
| 92 | #creates pts file for onshore DEM |
---|
| 93 | print "creates pts file for onshore DEM" |
---|
[5811] | 94 | dem2pts(onshore_dir_name ,use_cache=True,verbose=True) |
---|
[5810] | 95 | dem2pts(offshore_dir_name1 ,use_cache=True,verbose=True) |
---|
| 96 | dem2pts(offshore_dir_name2 ,use_cache=True,verbose=True) |
---|
| 97 | dem2pts(offshore_dir_name3 ,use_cache=True,verbose=True) |
---|
| 98 | dem2pts(offshore_dir_name4 ,use_cache=True,verbose=True) |
---|
| 99 | dem2pts(offshore_dir_name5 ,use_cache=True,verbose=True) |
---|
| 100 | dem2pts(offshore_dir_name6 ,use_cache=True,verbose=True) |
---|
| 101 | |
---|
| 102 | ###creates pts file for island DEM |
---|
| 103 | ##dem2pts(island_dir_name, use_cache=True, verbose=True) |
---|
| 104 | ##dem2pts(island_dir_name1, use_cache=True, verbose=True) |
---|
| 105 | ##dem2pts(island_dir_name2, use_cache=True, verbose=True) |
---|
| 106 | ##dem2pts(island_dir_name3, use_cache=True, verbose=True) |
---|
| 107 | ##dem2pts(island_dir_name4, use_cache=True, verbose=True) |
---|
| 108 | ##dem2pts(island_dir_name5, use_cache=True, verbose=True) |
---|
| 109 | ##dem2pts(island_dir_name6, use_cache=True, verbose=True) |
---|
| 110 | |
---|
| 111 | #------------------------------------------------------------------------------- |
---|
| 112 | # Combine datasets into project.combined_dir_name |
---|
| 113 | #------------------------------------------------------------------------------- |
---|
| 114 | |
---|
| 115 | print'create Geospatial data1 objects from topographies' |
---|
| 116 | G1 = Geospatial_data(file_name = onshore_dir_name + '.pts') |
---|
| 117 | print'create Geospatial data6 objects from topographies' |
---|
| 118 | G_coast = Geospatial_data(file_name = coast_in_dir_name) |
---|
| 119 | print'create Geospatial data7 objects from topographies' |
---|
| 120 | G_off = Geospatial_data(file_name = offshore_in_dir_name) |
---|
| 121 | print'create Geospatial data8 objects from topographies' |
---|
[5812] | 122 | G_off1 = Geospatial_data(file_name = offshore_dir_name1 + '.pts') |
---|
[5810] | 123 | print'create Geospatial data9 objects from topographies' |
---|
[5812] | 124 | G_off2 = Geospatial_data(file_name = offshore_dir_name2 + '.pts') |
---|
[5810] | 125 | print'create Geospatial data10 objects from topographies' |
---|
[5813] | 126 | G_off3 = Geospatial_data(file_name = offshore_dir_name3 + '.pts') |
---|
[5810] | 127 | print'create Geospatial data10 objects from topographies' |
---|
[5813] | 128 | G_off4 = Geospatial_data(file_name = offshore_dir_name4 + '.pts') |
---|
[5810] | 129 | print'create Geospatial data10 objects from topographies' |
---|
[5813] | 130 | G_off5 = Geospatial_data(file_name = offshore_dir_name5 + '.pts') |
---|
[5810] | 131 | print'create Geospatial data10 objects from topographies' |
---|
[5813] | 132 | G_off6 = Geospatial_data(file_name = offshore_dir_name6 + '.pts') |
---|
[5810] | 133 | print'add all geospatial objects' |
---|
[5811] | 134 | G = G1 + G_coast + G_off + G_off1 + G_off2 + G_off3 + G_off4 + G_off5 + G_off6 |
---|
[5810] | 135 | print'clip combined geospatial object by bounding polygon' |
---|
| 136 | G_clipped = G.clip(project.poly_all) |
---|
| 137 | |
---|
| 138 | print'export combined DEM file' |
---|
| 139 | if access(project.topographies_dir,F_OK) == 0: |
---|
| 140 | mkdir (project.topographies_dir) |
---|
| 141 | G_clipped.export_points_file(project.combined_dir_name + '.pts') |
---|
| 142 | #G_clipped.export_points_file(project.combined_dir_name + '.txt') #Use for comparision in ARC |
---|
| 143 | |
---|
| 144 | import sys |
---|
| 145 | sys.exit() |
---|