source: anuga_work/production/busselton/build_busselton_250m.py @ 5786

Last change on this file since 5786 was 5786, checked in by kristy, 16 years ago

Updated all scripts to reflect Perth format. Also added 250m scripts

File size: 3.2 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_250m.py
6Output: pts file stored in project_250m.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_250m.output_build_time_dir,__file__, 
38               dirname(project_250m.__file__)+sep+ project_250m.__name__+'.py' )
39
40start_screen_catcher(project_250m.output_build_time_dir)
41
42print 'USER:    ', project_250m.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#-------------------------------------------------------------------------------
51print "project_250m.poly_all", project_250m.poly_all
52print "project_250m.combined_dir_name", project_250m.combined_dir_name
53
54# input elevation directory filenames
55onshore_in_dir_name = project_250m.onshore_in_dir_name
56
57# output elevation directory filenames
58onshore_dir_name = project_250m.onshore_dir_name
59
60# creates DEM from asc data
61print "creates DEMs from asc data"
62convert_dem_from_ascii2netcdf(onshore_in_dir_name, basename_out=onshore_dir_name, use_cache=True, verbose=True)
63
64# creates pts file for onshore DEM
65print "creates pts file for onshore DEM"
66dem2pts(onshore_dir_name ,use_cache=True,verbose=True)
67
68# create onshore pts files
69print'create Geospatial data1 objects from topographies'
70G = Geospatial_data(file_name = onshore_dir_name + '.pts')
71
72#-------------------------------------------------------------------------------
73# Combine, clip and export dataset
74#-------------------------------------------------------------------------------
75
76print'clip combined geospatial object by bounding polygon'
77G_clipped = G.clip(project_250m.poly_all)
78
79print'export combined DEM file'
80if access(project_250m.topographies_dir,F_OK) == 0:
81    mkdir (project_250m.topographies_dir)
82G_clipped.export_points_file(project_250m.combined_dir_name + '.pts')
83#G_clipped.export_points_file(project_250m.combined_dir_name + '.txt') #Use for comparision in ARC
84
Note: See TracBrowser for help on using the repository browser.