[6286] | 1 | """Build the elevation data to run a tsunami inundation scenario |
---|
| 2 | for Busselton, WA, Australia. |
---|
| 3 | |
---|
| 4 | Input: elevation data from project.py |
---|
| 5 | Output: pts file stored in project.topographies_dir |
---|
| 6 | The run_busselton.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 | |
---|
| 17 | # Related major packages |
---|
| 18 | from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf |
---|
| 19 | from anuga.shallow_water.data_manager import dem2pts |
---|
| 20 | from anuga.shallow_water.data_manager import start_screen_catcher |
---|
| 21 | from anuga.shallow_water.data_manager import copy_code_files |
---|
| 22 | from anuga.geospatial_data.geospatial_data import Geospatial_data |
---|
| 23 | |
---|
| 24 | # Application specific imports |
---|
| 25 | import project # Definition of file names and polygons |
---|
| 26 | |
---|
| 27 | |
---|
| 28 | #------------------------------------------------------------------------------ |
---|
| 29 | # Copy scripts to time stamped output directory and capture screen |
---|
| 30 | # output to file |
---|
| 31 | #------------------------------------------------------------------------------ |
---|
| 32 | copy_code_files(project.output_build,__file__, |
---|
| 33 | os.path.dirname(project.__file__)+os.sep+\ |
---|
| 34 | project.__name__+'.py' ) |
---|
| 35 | start_screen_catcher(project.output_build) |
---|
| 36 | |
---|
| 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 | ##for filename in project.ascii_grid_filenames: |
---|
| 52 | ## absolute_filename = project.topographies_folder + os.sep + filename |
---|
| 53 | ## convert_dem_from_ascii2netcdf(absolute_filename, |
---|
| 54 | ## basename_out=absolute_filename, |
---|
| 55 | ## use_cache=True, |
---|
| 56 | ## verbose=True) |
---|
| 57 | ## dem2pts(absolute_filename, use_cache=True, verbose=True) |
---|
| 58 | ## |
---|
| 59 | ## geospatial_data[filename] = Geospatial_data(file_name=absolute_filename+'.pts', |
---|
| 60 | ## verbose=True) |
---|
| 61 | |
---|
| 62 | # Create Geospatial data from TXT files |
---|
| 63 | for filename in project.point_filenames: |
---|
| 64 | print filename |
---|
| 65 | absolute_filename = project.topographies_folder + os.sep + str(filename) |
---|
| 66 | geospatial_data[filename] = Geospatial_data(file_name=absolute_filename, |
---|
| 67 | verbose=True) |
---|
| 68 | |
---|
| 69 | |
---|
| 70 | #------------------------------------------------------------------------------- |
---|
| 71 | # Combine, clip and export dataset |
---|
| 72 | #------------------------------------------------------------------------------- |
---|
| 73 | |
---|
| 74 | print 'Add geospatial objects' # except', project.offshore_name5 |
---|
| 75 | G = None |
---|
| 76 | for key in geospatial_data: |
---|
| 77 | #if key == project.point_filenames[4] or key == project.ascii_filenames[1]: |
---|
| 78 | # # Skip these files for now |
---|
| 79 | # continue |
---|
| 80 | |
---|
| 81 | G += geospatial_data[key] |
---|
| 82 | |
---|
| 83 | print 'Clip combined geospatial data' |
---|
| 84 | ##G_clip = G.clip_outside(project.poly_aoi1) |
---|
| 85 | ##G_all = G_clip + geospatial_data[project.point_filenames[4]] |
---|
| 86 | #G_clipped = G_all.clip(project.poly_all) |
---|
| 87 | G_clip = G.clip(project.bounding_polygon) |
---|
| 88 | |
---|
| 89 | |
---|
| 90 | print 'Export combined DEM file' |
---|
[6288] | 91 | G_clip.export_points_file(project.combined_elevation + '.pts') |
---|
[6286] | 92 | print 'Do txt version too' |
---|
| 93 | # Use for comparision in ARC |
---|
[6288] | 94 | G_clip.export_points_file(project.combined_elevation + '.txt') |
---|
[6286] | 95 | |
---|