source: anuga_work/production/australia_ph2/ceduna/build_elevation.py @ 6359

Last change on this file since 6359 was 6359, checked in by myall, 16 years ago

phase 2 scripts for ceduna, eucla and esperance still not working...

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