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

Last change on this file since 6254 was 6254, checked in by ole, 15 years ago

Established work area for the development of a modelling demo

File size: 6.3 KB
RevLine 
[6253]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 import sep
17from os.path import dirname, basename
18from os import mkdir, access, F_OK
19from shutil import copy
20import time
21import sys
22
23# Related major packages
24from anuga.shallow_water import Domain
25from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf, dem2pts
26from anuga.geospatial_data.geospatial_data import *
27from anuga.shallow_water.data_manager import start_screen_catcher, copy_code_files
28
29# Application specific imports
30import project   # Definition of file names and polygons
31
32#------------------------------------------------------------------------------
33# Copy scripts to time stamped output directory and capture screen
34# output to file
35#------------------------------------------------------------------------------
36
37copy_code_files(project.output_build_time_dir,__file__, 
38               dirname(project.__file__)+sep+ project.__name__+'.py' )
39
40start_screen_catcher(project.output_build_time_dir)
41
42print 'USER:    ', project.user
43
44#-------------------------------------------------------------------------------
45# Preparation of topographic data
46#
47# Convert ASC 2 DEM 2 PTS using source data and store result in source data
48# Do for coarse and fine data
49# Fine pts file to be clipped to area of interest
50#-------------------------------------------------------------------------------
[6254]51
52# FIXME(Ole): Here's the Shark Bay way.
53
54## FIXME (Ole): Clip data by interior regions if possible
55#
56#print 'creating geospatial data objects from asc data (via dem and pts formats)'
57#for filename in project.ascii_grid_filenames:
58#    convert_dem_from_ascii2netcdf(filename,
59#                                  basename_out=filename,
60#                                  use_cache=True, verbose=True)
61#    dem2pts(filename, use_cache=True, verbose=True)#
62#
63#    geospatial_data += Geospatial_data(file_name=filename + '.pts',
64#                                       verbose=True)
65##
66#
67#print 'creating geospatial data objects from txt data'
68#for filename in project.point_filenames:
69#    geospatial_data += Geospatial_data(file_name=filename + '.txt',
70#                                       verbose=True)
71#
72#
73#print 'clip combined geospatial object by bounding polygon'
74#G = geospatial_data.clip(project.bounding_polygon)
75#
76#
77#print 'export combined geospatial data'
78#if access(project.topographies_dir, F_OK) == 0:
79#    mkdir (project.topographies_dir)
80#G.export_points_file(project.combined_dir_name + '.pts')
81
82
83
[6253]84print "project.poly_all", project.poly_all
85print "project.combined_dir_name", project.combined_dir_name
86
87# input elevation directory filenames
88onshore_in_dir_name = project.onshore_in_dir_name
89coast_in_dir_name = project.coast_in_dir_name
90coast_in_dir_name1 = project.coast_in_dir_name1
91offshore_in_dir_name = project.offshore_in_dir_name
92offshore_in_dir_name1 = project.offshore_in_dir_name1
93offshore_in_dir_name2 = project.offshore_in_dir_name2
94offshore_in_dir_name3 = project.offshore_in_dir_name3
95offshore_in_dir_name4 = project.offshore_in_dir_name4
96offshore_in_dir_name5 = project.offshore_in_dir_name5
97
98# output elevation directory filenames
99onshore_dir_name = project.onshore_dir_name
100coast_dir_name = project.coast_dir_name
101coast_dir_name1 = project.coast_dir_name1
102offshore_dir_name = project.offshore_dir_name
103offshore_dir_name1 = project.offshore_dir_name1
104offshore_dir_name2 = project.offshore_dir_name2
105offshore_dir_name3 = project.offshore_dir_name3
106offshore_dir_name4 = project.offshore_dir_name4
107offshore_dir_name5 = project.offshore_dir_name5
108# creates DEM from asc data
109print "creates DEMs from asc data"
110convert_dem_from_ascii2netcdf(onshore_in_dir_name, basename_out=onshore_dir_name, use_cache=True, verbose=True)
111convert_dem_from_ascii2netcdf(offshore_in_dir_name5, basename_out=offshore_dir_name5, use_cache=True, verbose=True)
112
113# creates pts file for onshore DEM
114print "creates pts file for onshore DEM"
115dem2pts(onshore_dir_name ,use_cache=True,verbose=True)
116dem2pts(offshore_dir_name5 ,use_cache=True,verbose=True)
117
118# create onshore pts files
119print'create Geospatial data1 objects from topographies'
120G1 = Geospatial_data(file_name = onshore_dir_name + '.pts')
121
122# create coastal and offshore pts files
123print'create Geospatial data2 objects from topographies'
124G2 = Geospatial_data(file_name = coast_in_dir_name)
125print'create Geospatial data3 objects from topographies'
126G3 = Geospatial_data(file_name = coast_in_dir_name1)
127print'create Geospatial data4 objects from topographies'
128G_off = Geospatial_data(file_name = offshore_in_dir_name)
129print'create Geospatial data5 objects from topographies'
130G_off1 = Geospatial_data(file_name = offshore_in_dir_name1)
131print'create Geospatial data6 objects from topographies'
132G_off2 = Geospatial_data(file_name = offshore_in_dir_name2)
133print'create Geospatial data7 objects from topographies'
134G_off3 = Geospatial_data(file_name = offshore_in_dir_name3)
135print'create Geospatial data8 objects from topographies'
136G_off4 = Geospatial_data(file_name = offshore_in_dir_name4)
137print'create Geospatial data9 objects from topographies'
138G_off5 = Geospatial_data(file_name = offshore_dir_name5 + '.pts')
139
140
141#-------------------------------------------------------------------------------
142# Combine, clip and export dataset
143#-------------------------------------------------------------------------------
144
145print'add all geospatial objects'
146G = G1 + G2 + G3 + G_off + G_off1 + G_off2 + G_off3 + G_off4
147
148print'clip combined geospatial object by bounding polygon'
149G_clip = G.clip_outside(project.poly_aoi1)
150G_all = G_clip + G_off5
151G_clipped = G_all.clip(project.poly_all)
152
153print'export combined DEM file'
154if access(project.topographies_dir,F_OK) == 0:
155    mkdir (project.topographies_dir)
156G_clipped.export_points_file(project.combined_dir_name + '.pts')
157G_clipped.export_points_file(project.combined_dir_name + '.txt') #Use for comparision in ARC
158
Note: See TracBrowser for help on using the repository browser.