source: anuga_work/production/bunbury_storm_surge_2009/build_elevation.py @ 7658

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

updates to bunbury storm surge model

File size: 6.5 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 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
84keylist = project.ascii_grid_filenames + project.point_filenames
85
86G1 = geospatial_data[keylist[0]].clip_outside(project.elevation_clip_box)
87G2 = geospatial_data[keylist[1]].clip(project.elevation_clip_box)
88G3 = geospatial_data[keylist[2]]
89G4 = geospatial_data[keylist[3]]
90G5 = geospatial_data[keylist[4]]
91G6 = geospatial_data[keylist[5]]
92G7 = geospatial_data[keylist[6]]
93G8 = geospatial_data[keylist[7]]
94G9 = geospatial_data[keylist[8]]
95G10 = geospatial_data[keylist[9]]
96#G11 = geospatial_data[keylist[10]]
97
98
99# G = None
100
101# for i, key in enumerate(keylist):
102    # if key==keylist[0]:
103        # G[i] = geospatial_data[key].clip_outside(project.elevation_clip_box)
104    # elif key==keylist[1]:
105        # G[i] = geospatial_data[key].clip(project.elevation_clip_box)
106    # else
107        # G+ = geospatial_data[key]
108       
109#~ G_list = [G1, G2, G3, G4, G5, G6, G7, G8, G9, G10]
110
111#~ print 'G1', G1.attributes
112#~ print 'G2', G2.attributes
113#~ print 'G3', G3.attributes
114#~ print 'G4', G4.attributes
115#~ print 'G5', G5.attributes
116#~ print 'G6', G6.attributes
117#~ print 'G7', G7.attributes
118#~ print 'G8', G8.attributes
119#~ print 'G9', G9.attributes
120#~ print 'G10', G10.attributes
121
122G = G1 + G2 + G3 + G4 + G5 + G6 + G7 + G8 + G9 + G10 #+ G11
123
124# G = None
125# for key in geospatial_data:
126    # if key == project.ascii_grid_filenames[0]:
127                # G == geospatial_data[key].clip_outside(project.elevation_clip_box)
128        # elif key == project.ascii_grid_filenames[1]:
129                # G = geospatial_data[key].clip(project.elevation_clip_box)
130    # continue
131   
132    # G += geospatial_data[key]
133
134print 'Export combined DEM file'
135G.export_points_file(project.combined_elevation + '.pts')
136
137print 'Do txt version too'
138# Use for comparision in ARC
139#~ for g in G_list:
140        #~ for i in keylist:
141        #~ g.export_points_file(join(project.topographies_folder, (str(keylist[i]) +'_export.txt')))
142
143G.export_points_file(project.combined_elevation + '.txt')
144
145print 'Export individual text files'
146#for i in xrange(len(keylist)):
147#    G[i+1].export_points_file(join(project.topographies_folder, keylist[i] +'_export.txt'))
148 
149G1.export_points_file(join(project.topographies_folder, keylist[0] +'_export.txt'))
150G2.export_points_file(join(project.topographies_folder, keylist[1] +'_export.txt'))
151G3.export_points_file(join(project.topographies_folder, keylist[2] +'_export.txt'))
152G4.export_points_file(join(project.topographies_folder, keylist[3] +'_export.txt'))
153G5.export_points_file(join(project.topographies_folder, keylist[4] +'_export.txt'))
154G6.export_points_file(join(project.topographies_folder, keylist[5] +'_export.txt'))
155G7.export_points_file(join(project.topographies_folder, keylist[6] +'_export.txt'))
156G8.export_points_file(join(project.topographies_folder, keylist[7] +'_export.txt'))
157G9.export_points_file(join(project.topographies_folder, keylist[8] +'_export.txt'))
158G10.export_points_file(join(project.topographies_folder, keylist[9] +'_export.txt'))
159#G11.export_points_file(join(project.topographies_folder, keylist[10] +'_export.txt'))
Note: See TracBrowser for help on using the repository browser.