source: anuga_work/production/wa/geraldton_2009/build_elevation.py @ 7855

Last change on this file since 7855 was 6569, checked in by kristy, 16 years ago

new scripts that encompass all changes to this date

File size: 4.4 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
24from anuga.utilities.polygon import read_polygon
25
26# Application specific imports
27from setup_model import project   # Definition of file names and polygons
28
29
30#------------------------------------------------------------------------------
31# Copy scripts to time stamped output directory and capture screen
32# output to file
33#------------------------------------------------------------------------------
34copy_code_files(project.output_build,__file__, 
35                os.path.dirname(project.__file__)+os.sep+\
36                project.__name__+'.py' )
37start_screen_catcher(project.output_build)
38
39
40#------------------------------------------------------------------------------
41# Preparation of topographic data
42#
43# Convert ASC 2 DEM 2 PTS using source data and store result in source data
44# Do for coarse and fine data
45# Fine pts file to be clipped to area of interest
46#------------------------------------------------------------------------------
47
48print 'project.bounding_polygon', project.bounding_polygon
49print 'project.combined_elevation_basename', project.combined_elevation_basename
50
51# Create Geospatial data from ASCII files
52geospatial_data = {}
53if not project.ascii_grid_filenames == []:
54    for filename in project.ascii_grid_filenames:
55        absolute_filename = join(project.topographies_folder, filename)
56        convert_dem_from_ascii2netcdf(absolute_filename,
57                                      basename_out=absolute_filename,
58                                      use_cache=True,
59                                      verbose=True)
60        dem2pts(absolute_filename, use_cache=True, verbose=True)
61
62        G_grid = Geospatial_data(file_name=absolute_filename+'.pts',
63                                                    verbose=True)
64        print 'Clip geospatial object'
65        geospatial_data[filename] = G_grid.clip(project.bounding_polygon)
66
67# Create Geospatial data from TXT files
68if not project.point_filenames == []:
69    for filename in project.point_filenames:
70        absolute_filename = join(project.topographies_folder, filename)
71        G_points = Geospatial_data(file_name=absolute_filename,
72                                                    verbose=True)
73        print 'Clip geospatial object'
74        geospatial_data[filename] = G_points.clip(project.bounding_polygon)
75
76#-------------------------------------------------------------------------------
77# Combine, clip and export dataset
78#-------------------------------------------------------------------------------
79
80extent_polygons = []
81for extent_polygon_filename in project.extent_polygon_filenames:
82    p = read_polygon(join(project.polygons_folder, extent_polygon_filename))
83    extent_polygons.append(p)
84   
85print 'Add geospatial objects except ', project.point_filenames[2]
86G = None
87for key in geospatial_data:
88    if key == project.point_filenames[2]:
89        D = geospatial_data[key]
90        D = D.clip(extent_polygons[0])
91        G += D
92    elif key == project.point_filenames[0] or key == project.point_filenames[1] or key == project.point_filenames[3] or key == project.point_filenames[4]:
93        D = geospatial_data[key]
94        D = D.clip(extent_polygons[1])
95        G += D
96    elif key == project.ascii_grid_filenames[0] or key == project.ascii_grid_filenames[1]:
97        D = geospatial_data[key]
98        D = D.clip(extent_polygons[2])
99        G += D
100
101
102print 'Export combined DEM file'
103G.export_points_file(project.combined_elevation + '.pts')
104print 'Do txt version too'
105# Use for comparision in ARC
106G.export_points_file(project.combined_elevation + '.txt')
107
Note: See TracBrowser for help on using the repository browser.