source: anuga_validation/automated_validation_tests/patong_beach_validation/build_elevation.py @ 7750

Last change on this file since 7750 was 6844, checked in by rwilson, 15 years ago

Add Patong beach validation in main trunk.

  • Property svn:executable set to *
File size: 4.1 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 = {}
53
54for 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
68
69for 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#-------------------------------------------------------------------------------
79extent_polygons = []
80for extent_polygon_filename in project.extent_polygon_filenames:
81    p = read_polygon(join(project.polygons_folder, extent_polygon_filename))
82    extent_polygons.append(p)
83   
84print 'Add geospatial objects' 
85G = None
86for key in geospatial_data:
87    if key == project.point_filenames[0] or key == project.point_filenames[1]:
88        G += geospatial_data[key]
89    elif key == project.point_filenames[2]:
90        D = geospatial_data[key]
91        D = D.clip_outside(extent_polygons[0])
92        D = D.clip_outside(extent_polygons[1])
93        G += D
94    elif key == project.point_filenames[3]:
95        D = geospatial_data[key]
96        D = D.clip_outside(extent_polygons[2])
97        G += D
98
99print 'Export combined DEM file'
100G.export_points_file(project.combined_elevation + '.pts')
101print 'Do txt version too'
102# Use for comparision in ARC
103G.export_points_file(project.combined_elevation + '.txt')
104
Note: See TracBrowser for help on using the repository browser.