source: anuga_work/production/wa/mandurah_2009/build_elevation.py @ 7557

Last change on this file since 7557 was 7557, checked in by monica, 14 years ago

updated data

File size: 3.9 KB
Line 
1"""Build the elevation data to run a tsunami inundation scenario
2for Busselton, WA, 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#------------------------------------------------------------------------------
40# Preparation of topographic data
41#
42# Convert ASC 2 DEM 2 PTS using source data and store result in source data
43# Do for coarse and fine data
44# Fine pts file to be clipped to area of interest
45#------------------------------------------------------------------------------
46
47print 'project.bounding_polygon', project.bounding_polygon
48print 'project.combined_elevation_basename', project.combined_elevation_basename
49
50# Create Geospatial data from ASCII files
51geospatial_data = {}
52if not project.ascii_grid_filenames == []:
53    for filename in project.ascii_grid_filenames:
54        absolute_filename = join(project.topographies_folder, filename)
55        convert_dem_from_ascii2netcdf(absolute_filename,
56                                      basename_out=absolute_filename,
57                                      use_cache=True,
58                                      verbose=True)
59        dem2pts(absolute_filename, use_cache=True, verbose=True)
60
61        G_grid = Geospatial_data(file_name=absolute_filename+'.pts',
62                                                    verbose=True)
63        print 'Clip geospatial object'
64        geospatial_data[filename] = G_grid.clip(project.bounding_polygon)
65
66#Create Geospatial data from TXT files
67if not project.point_filenames == []:
68    for filename in project.point_filenames:
69        absolute_filename = join(project.topographies_folder, filename)
70        G_points = Geospatial_data(file_name=absolute_filename,
71                                                    verbose=True)
72        print 'Clip geospatial object'
73        geospatial_data[filename] = G_points.clip(project.bounding_polygon)
74
75#-------------------------------------------------------------------------------
76# Combine, clip and export dataset
77#-------------------------------------------------------------------------------
78
79print 'Add geospatial objects' # except', project.offshore_name5
80G = None
81G1 = None
82G2 = None
83for key in geospatial_data:
84    if key == 'extract_m_peel':
85        G1 = geospatial_data[key].clip(project.clip_box)
86    if key == 'm_peel_10m':
87        G2 = geospatial_data[key].clip_outside(project.clip_box)
88
89    G += geospatial_data[key]
90G = G + G1 + G2
91   
92print 'Export combined DEM file'
93G.export_points_file(project.combined_elevation + '.pts')
94print 'Do txt version too'
95# Use for comparision in ARC
96G.export_points_file(project.combined_elevation + '.txt')
Note: See TracBrowser for help on using the repository browser.