source: anuga_work/production/mandurah_storm_surge_2009/gems_comparison/build_elevation_GEMS.py @ 7647

Last change on this file since 7647 was 7647, checked in by fountain, 14 years ago

updates to mandurah and bunbury storm surge models

File size: 4.2 KB
Line 
1"""Build the elevation data to run a storm surge inundation scenario
2for Mandurah, Western Australia, Australia.
3
4Input: elevation data from project.py
5Output: pts file stored in project.topographies_dir
6The run_model.py is reliant on the output of this script.
7
8"""
9
10#------------------------------------------------------------------------------
11# Import necessary modules
12#------------------------------------------------------------------------------
13
14# Standard modules
15import os
16from os.path import join
17
18# Related major packages
19from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf
20from anuga.shallow_water.data_manager import dem2pts
21from anuga.shallow_water.data_manager import start_screen_catcher
22from anuga.shallow_water.data_manager import copy_code_files
23from anuga.geospatial_data.geospatial_data import Geospatial_data
24
25# Application specific imports
26from setup_model_GEMS import project_GEMS   # 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#------------------------------------------------------------------------------
33copy_code_files(project_GEMS.output_build,__file__, 
34                                os.path.dirname(project_GEMS.__file__)+os.sep+\
35                                project_GEMS.__name__+'.py' )
36start_screen_catcher(project_GEMS.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
46print 'project.bounding_polygon', project_GEMS.bounding_polygon
47print 'project.combined_elevation_basename', project_GEMS.combined_elevation_basename
48
49# Create Geospatial data from ASCII files
50geospatial_data = {}
51if not project_GEMS.ascii_grid_filenames == []:
52    for filename in project_GEMS.ascii_grid_filenames:
53        absolute_filename = join(project_GEMS.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_GEMS.bounding_polygon)
64
65# Create Geospatial data from TXT files
66if not project_GEMS.point_filenames == []:
67    for filename in project_GEMS.point_filenames:
68        absolute_filename = join(project_GEMS.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_GEMS.bounding_polygon)
73
74print 'Geospatial data objects: ', geospatial_data.keys()
75
76#-------------------------------------------------------------------------------
77# Combine, clip and export dataset
78#-------------------------------------------------------------------------------
79
80print 'Add geospatial objects' # except', project_GEMS.offshore_name5
81
82print 'project_GEMS.elevation_clip_box', project_GEMS.elevation_clip_box
83
84G1 = geospatial_data['man_25m_prj'].clip(project_GEMS.elevation_clip_box)
85G2 = geospatial_data['man10m_ss'].clip_outside(project_GEMS.elevation_clip_box)
86
87G = G1 + G2
88       
89print 'Export combined DEM file'
90G.export_points_file(project_GEMS.combined_elevation + '.pts')
91
92print 'Do txt version too'
93# Use for comparision in ARC
94G.export_points_file(project_GEMS.combined_elevation + '.txt')
95G1.export_points_file(join(project_GEMS.topographies_folder, (str(project_GEMS.ascii_grid_filenames[0]) +'_export.txt')))
96G2.export_points_file(join(project_GEMS.topographies_folder, (str(project_GEMS.ascii_grid_filenames[1]) + '_export.txt')))
Note: See TracBrowser for help on using the repository browser.