source: anuga_work/production/mandurah_storm_surge_2009/build_elevation.py @ 7627

Last change on this file since 7627 was 7627, checked in by fountain, 13 years ago

updates to mandurah storm surge scenario

File size: 4.4 KB
RevLine 
[7614]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 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#------------------------------------------------------------------------------
33copy_code_files(project.output_build,__file__, 
34                                os.path.dirname(project.__file__)+os.sep+\
35                                project.__name__+'.py' )
36start_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
46print 'project.bounding_polygon', project.bounding_polygon
47print 'project.combined_elevation_basename', project.combined_elevation_basename
48
49# Create Geospatial data from ASCII files
50geospatial_data = {}
51if 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
66if 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
74print 'Geospatial data objects: ', geospatial_data.keys()
75
76#-------------------------------------------------------------------------------
77# Combine, clip and export dataset
78#-------------------------------------------------------------------------------
79
80print 'Add geospatial objects' # except', project.offshore_name5
81
82print 'project.elevation_clip_box', project.elevation_clip_box
83
[7627]84# G1 = geospatial_data['Man_25m'].clip(project.elevation_clip_box)
85# G2 = geospatial_data['ph10m_ss'].clip_outside(project.elevation_clip_box)
[7614]86
[7627]87G1 = geospatial_data['m_peel_aoi'].clip(project.elevation_clip_box)
88G2 = geospatial_data['ph10m_ss'].clip_outside(project.elevation_clip_box)
89G3 = geospatial_data['MA-46893-SNDS_AHD.csv']
90G4 = geospatial_data['MS0205HY_AHD.csv']
91G5 = geospatial_data['MS0404_AHD.csv']
92G6 = geospatial_data['YU0403HY_AHD.csv']
93G7 = geospatial_data['original_data_ss.csv']
94
95G = G1 + G2 + G3 + G4 + G5 + G6 + G7
[7614]96       
97print 'Export combined DEM file'
98G.export_points_file(project.combined_elevation + '.pts')
99
100print 'Do txt version too'
101# Use for comparision in ARC
102G.export_points_file(project.combined_elevation + '.txt')
103G1.export_points_file(join(project.topographies_folder, (str(project.ascii_grid_filenames[0]) +'_export.txt')))
104G2.export_points_file(join(project.topographies_folder, (str(project.ascii_grid_filenames[1]) + '_export.txt')))
Note: See TracBrowser for help on using the repository browser.