source: anuga_work/production/australia_ph2/albany/build_elevation.py @ 6484

Last change on this file since 6484 was 6484, checked in by myall, 14 years ago

scripts for running albany model

File size: 3.8 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        geospatial_data[filename] = Geospatial_data(file_name=absolute_filename+'.pts',
62                                                    verbose=True)
63
64# Create Geospatial data from TXT files
65if not project.point_filenames == []:
66    for filename in project.point_filenames:
67        absolute_filename = join(project.topographies_folder, filename)
68        geospatial_data[filename] = Geospatial_data(file_name=absolute_filename,
69                                                    verbose=True)
70
71#-------------------------------------------------------------------------------
72# Combine, clip and export dataset
73#-------------------------------------------------------------------------------
74
75print 'Add geospatial objects' # except', project.offshore_name5
76G = None
77for key in geospatial_data:
78    #if key == project.point_filenames[4] or key == project.ascii_filenames[1]:
79    #    # Skip these files for now
80    #    continue
81   
82    G += geospatial_data[key]
83       
84print 'Clip combined geospatial data'
85##G_clip = G.clip_outside(project.poly_aoi1)
86##G_all = G_clip + geospatial_data[project.point_filenames[4]]
87#G_clipped = G_all.clip(project.poly_all)
88G_clip = G.clip(project.bounding_polygon)
89
90
91print 'Export combined DEM file'
92G_clip.export_points_file(project.combined_elevation + '.pts')
93print 'Do txt version too'
94# Use for comparision in ARC
95G_clip.export_points_file(project.combined_elevation + '.txt')
96
Note: See TracBrowser for help on using the repository browser.