source: anuga_work/production/perth/build_perth.py @ 6748

Last change on this file since 6748 was 5784, checked in by sexton, 16 years ago

updated commments for perth scripts

File size: 5.9 KB
RevLine 
[5761]1"""
[5784]2Script for building the elevation data to run a tsunami inundation scenario
[5761]3for Perth, WA, Australia.
[4077]4
[5761]5Input: elevation data from project.py
6Output: pts file stored in project.topographies_dir
[5759]7The run_perth.py is reliant on the output of this script.
8
[4077]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 *
[5360]27from anuga.shallow_water.data_manager import start_screen_catcher, copy_code_files
[4077]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#-------------------------------------------------------------------------------
[5784]51print "project.poly_all", project.poly_all
52print "project.combined_dir_name", project.combined_dir_name
[4091]53
[5759]54# input elevation directory filenames
[4080]55onshore_in_dir_name = project.onshore_in_dir_name
56coast_in_dir_name = project.coast_in_dir_name
57island_in_dir_name = project.island_in_dir_name
58island_in_dir_name1 = project.island_in_dir_name1
59island_in_dir_name2 = project.island_in_dir_name2
60island_in_dir_name3 = project.island_in_dir_name3
61offshore_in_dir_name = project.offshore_in_dir_name
[5524]62offshore_in_dir_name1 = project.offshore_in_dir_name1
63offshore_in_dir_name2 = project.offshore_in_dir_name2
[5753]64offshore_in_dir_name3 = project.offshore_in_dir_name3
[4077]65
[5759]66# output elevation directory filenames
[4077]67onshore_dir_name = project.onshore_dir_name
68coast_dir_name = project.coast_dir_name
[4080]69island_dir_name = project.island_dir_name
70island_dir_name1 = project.island_dir_name1
71island_dir_name2 = project.island_dir_name2
72island_dir_name3 = project.island_dir_name3
[4077]73offshore_dir_name = project.offshore_dir_name
[5524]74offshore_dir_name1 = project.offshore_dir_name1
75offshore_dir_name2 = project.offshore_dir_name2
[5753]76offshore_dir_name3 = project.offshore_dir_name3
[4077]77
78# creates DEM from asc data
[4091]79print "creates DEMs from asc data"
[4080]80convert_dem_from_ascii2netcdf(onshore_in_dir_name, basename_out=onshore_dir_name, use_cache=True, verbose=True)
81convert_dem_from_ascii2netcdf(island_in_dir_name, basename_out=island_dir_name, use_cache=True, verbose=True)
82convert_dem_from_ascii2netcdf(island_in_dir_name1, basename_out=island_dir_name1, use_cache=True, verbose=True)
83convert_dem_from_ascii2netcdf(island_in_dir_name2, basename_out=island_dir_name2, use_cache=True, verbose=True)
84convert_dem_from_ascii2netcdf(island_in_dir_name3, basename_out=island_dir_name3, use_cache=True, verbose=True)
[4077]85
[5784]86# creates pts file for onshore DEM
[4091]87print "creates pts file for onshore DEM"
[5759]88dem2pts(onshore_dir_name ,use_cache=True,verbose=True)
[4077]89
[5784]90# creates pts file for island DEM
[4080]91dem2pts(island_dir_name, use_cache=True, verbose=True)
92dem2pts(island_dir_name1, use_cache=True, verbose=True)
93dem2pts(island_dir_name2, use_cache=True, verbose=True)
94dem2pts(island_dir_name3, use_cache=True, verbose=True)
[4077]95
[5784]96# create onshore pts files
[4091]97print'create Geospatial data1 objects from topographies'
[4080]98G1 = Geospatial_data(file_name = onshore_dir_name + '.pts')
[4091]99print'create Geospatial data2 objects from topographies'
[5759]100G2 = Geospatial_data(file_name = island_dir_name + '.pts')
[4091]101print'create Geospatial data3 objects from topographies'
[5759]102G3 = Geospatial_data(file_name = island_dir_name1 + '.pts')
[4091]103print'create Geospatial data4 objects from topographies'
[5759]104G4 = Geospatial_data(file_name = island_dir_name2 + '.pts')
[4091]105print'create Geospatial data5 objects from topographies'
[5759]106G5 = Geospatial_data(file_name = island_dir_name3 + '.pts')
[5784]107
108# create coastal and offshore pts files
[4091]109print'create Geospatial data6 objects from topographies'
[5759]110G_coast = Geospatial_data(file_name = coast_in_dir_name)
[4091]111print'create Geospatial data7 objects from topographies'
[5759]112G_off = Geospatial_data(file_name = offshore_in_dir_name)
[5753]113print'create Geospatial data8 objects from topographies'
[5759]114G_off1 = Geospatial_data(file_name = offshore_in_dir_name1)
[5753]115print'create Geospatial data9 objects from topographies'
[5759]116G_off2 = Geospatial_data(file_name = offshore_in_dir_name2)
[5753]117print'create Geospatial data10 objects from topographies'
[5759]118G_off3 = Geospatial_data(file_name = offshore_in_dir_name3)
[4077]119
[5784]120#-------------------------------------------------------------------------------
121# Combine, clip and export dataset
122#-------------------------------------------------------------------------------
[5360]123
[4077]124print'add all geospatial objects'
[5759]125G = G1 + G2 + G3 + G4 + G5 + G_coast + G_off + G_off1 + G_off2 + G_off3
[4077]126
[4080]127print'clip combined geospatial object by bounding polygon'
[5386]128G_clipped = G.clip(project.poly_all)
[4077]129
130print'export combined DEM file'
131if access(project.topographies_dir,F_OK) == 0:
132    mkdir (project.topographies_dir)
[5784]133G_clipped.export_points_file(project.combined_dir_name + '.pts')
134#G_clipped.export_points_file(project.combined_dir_name + '.txt') #Use for comparision in ARC
[4077]135
Note: See TracBrowser for help on using the repository browser.