source: anuga_work/production/busselton/standardised_version/build_busselton.py @ 6261

Last change on this file since 6261 was 6261, checked in by ole, 10 years ago

Looped elevation data conversions

File size: 6.5 KB
Line 
1"""
2Script for building the elevation data to run a tsunami inundation scenario
3for busselton, WA, Australia.
4
5Input: elevation data from project.py
6Output: pts file stored in project.topographies_dir
7The run_busselton.py is reliant on the output of this script.
8
9"""
10
11#------------------------------------------------------------------------------
12# Import necessary modules
13#------------------------------------------------------------------------------
14
15# Standard modules
16from os.path import dirname
17from os import mkdir, access, F_OK
18
19# Related major packages
20from anuga.shallow_water import Domain
21from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf, dem2pts
22from anuga.geospatial_data.geospatial_data import *
23from anuga.shallow_water.data_manager import start_screen_catcher, copy_code_files
24
25# Application specific imports
26import project   # Definition of file names and polygons
27
28#------------------------------------------------------------------------------
29# Copy scripts to time stamped output directory and capture screen
30# output to file
31#------------------------------------------------------------------------------
32
33#copy_code_files(project.output_build_time_dir,__file__,
34#               dirname(project.__file__)+os.sep+ project.__name__+'.py' )
35#start_screen_catcher(project.output_build_time_dir)
36
37
38#-------------------------------------------------------------------------------
39# Preparation of topographic data
40#
41# Convert ASC 2 DEM 2 PTS using source data and store result in source data
42# Do for coarse and fine data
43# Fine pts file to be clipped to area of interest
44#-------------------------------------------------------------------------------
45
46print 'project.poly_all', project.poly_all
47print 'project.combined_dir_name', project.combined_dir_name
48
49# input elevation directory filenames
50#onshore_in_dir_name = project.onshore_in_dir_name
51#coast_in_dir_name = project.coast_in_dir_name
52#coast_in_dir_name1 = project.coast_in_dir_name1
53#offshore_in_dir_name = project.offshore_in_dir_name
54#offshore_in_dir_name1 = project.offshore_in_dir_name1
55#offshore_in_dir_name2 = project.offshore_in_dir_name2
56#offshore_in_dir_name3 = project.offshore_in_dir_name3
57#offshore_in_dir_name4 = project.offshore_in_dir_name4
58#offshore_in_dir_name5 = project.offshore_in_dir_name5
59
60# output elevation directory filenames
61#onshore_dir_name = project.onshore_dir_name
62#coast_dir_name = project.coast_dir_name
63#coast_dir_name1 = project.coast_dir_name1
64#offshore_dir_name = project.offshore_dir_name
65#offshore_dir_name1 = project.offshore_dir_name1
66#offshore_dir_name2 = project.offshore_dir_name2
67#offshore_dir_name3 = project.offshore_dir_name3
68#offshore_dir_name4 = project.offshore_dir_name4
69#offshore_dir_name5 = project.offshore_dir_name5
70
71
72# Create Geospatial data from ASCII files
73geospatial_data = {}
74for filename in project.ascii_grid_filenames:
75    absolute_filename = project.topographies_in_dir + filename
76    convert_dem_from_ascii2netcdf(absolute_filename,
77                                  basename_out=absolute_filename,
78                                  use_cache=True,
79                                  verbose=True)
80    dem2pts(absolute_filename, use_cache=True, verbose=True)
81
82    geospatial_data[filename] = Geospatial_data(file_name=absolute_filename+'.pts',
83                                                verbose=True)
84
85# Create Geospatial data from TXT files
86for filename in project.point_filenames:
87    absolute_filename = project.topographies_in_dir + filename
88    geospatial_data[filename] = Geospatial_data(file_name=absolute_filename,
89                                                verbose=True)
90
91
92#print "creates DEMs from asc data"
93#print onshore_in_dir_name, onshore_dir_name
94#print offshore_in_dir_name5, offshore_dir_name5
95#
96#convert_dem_from_ascii2netcdf(onshore_in_dir_name, basename_out=onshore_dir_name, use_cache=True, verbose=True)
97#
98#convert_dem_from_ascii2netcdf(offshore_in_dir_name5, basename_out=offshore_dir_name5, use_cache=True, verbose=True)#
99##
100#
101# creates pts file for onshore DEM
102#print "creates pts file for onshore DEM"
103#dem2pts(onshore_dir_name ,use_cache=True,verbose=True)
104#dem2pts(offshore_dir_name5 ,use_cache=True,verbose=True)
105
106# create onshore pts files
107#print'create Geospatial data1 objects from topographies'
108#G1 = Geospatial_data(file_name = onshore_dir_name + '.pts')
109
110# create coastal and offshore pts files
111#print'create Geospatial data2 objects from topographies'
112#G2 = Geospatial_data(file_name = coast_in_dir_name)
113#print'create Geospatial data3 objects from topographies'
114#G3 = Geospatial_data(file_name = coast_in_dir_name1)
115#print'create Geospatial data4 objects from topographies'
116#G_off = Geospatial_data(file_name = offshore_in_dir_name)
117#print'create Geospatial data5 objects from topographies'
118#G_off1 = Geospatial_data(file_name = offshore_in_dir_name1)
119#print'create Geospatial data6 objects from topographies'
120#G_off2 = Geospatial_data(file_name = offshore_in_dir_name2)
121#print'create Geospatial data7 objects from topographies'
122#G_off3 = Geospatial_data(file_name = offshore_in_dir_name3)
123#print'create Geospatial data8 objects from topographies'
124#G_off4 = Geospatial_data(file_name = offshore_in_dir_name4)
125#print'create Geospatial data9 objects from topographies'
126#G_off5 = Geospatial_data(file_name = offshore_dir_name5 + '.pts')
127
128
129#-------------------------------------------------------------------------------
130# Combine, clip and export dataset
131#-------------------------------------------------------------------------------
132
133#print'add all geospatial objects'
134#G = G1 + G2 + G3 + G_off + G_off1 + G_off2 + G_off3 + G_off4
135
136print 'Add geospatial objects except', project.offshore_name5
137print geospatial_data.keys()
138G = None
139for key in geospatial_data:
140    if key != project.offshore_name5:
141        G += geospatial_data[key]
142       
143print 'Clip (outside) combined geospatial object except %s by area of interest.'
144G_clip = G.clip_outside(project.poly_aoi1)
145
146print 'Clip combined geospatial object by bounding polygon'
147G_all = G_clip + geospatial_data[project.offshore_name5]   # G_off5
148G_clipped = G_all.clip(project.poly_all)
149
150print 'Export combined DEM file'
151#if access(project.topographies_dir,F_OK) == 0:
152#    mkdir (project.topographies_dir)
153G_clipped.export_points_file(project.combined_dir_name + '.pts')
154print 'Do txt version too'
155G_clipped.export_points_file(project.combined_dir_name + '.txt') # Use for comparision in ARC
156
Note: See TracBrowser for help on using the repository browser.