source: anuga_work/production/patong/build_patong.py @ 7613

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

Attempt to get Patong running with wider extent to minimise boundary effects.
Mesh is still as original though. This needs to fixed.

File size: 3.8 KB
Line 
1"""
2Script for building the elevation data to run a tsunami inundation scenario
3for patong, Thailand.
4
5Input: elevation data from project.py
6Output: pts file stored in project.topographies_dir
7The run_patong.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#-------------------------------------------------------------------------------
51print "project.poly_all", project.poly_all
52print "project.combined_dir_name", project.combined_dir_name
53
54# input elevation directory filenames
55elevation_in_dir_name1 = project.elevation_in_dir_name1
56elevation_in_dir_name2 = project.elevation_in_dir_name2
57elevation_in_dir_name3 = project.elevation_in_dir_name3
58elevation_in_dir_name4 = project.elevation_in_dir_name4
59
60# output elevation directory filenames
61elevation_dir_name1 = project.elevation_dir_name1
62elevation_dir_name2 = project.elevation_dir_name2
63elevation_dir_name3 = project.elevation_dir_name3
64elevation_dir_name4 = project.elevation_dir_name4
65
66# create elevation pts files
67print'create Geospatial data1 objects from topographies'
68G1 = Geospatial_data(file_name = elevation_in_dir_name1)
69print'create Geospatial data2 objects from topographies'
70G2 = Geospatial_data(file_name = elevation_in_dir_name2)
71print'create Geospatial data3 objects from topographies'
72G3 = Geospatial_data(file_name = elevation_in_dir_name3)
73print'create Geospatial data4 objects from topographies'
74G4 = Geospatial_data(file_name = elevation_in_dir_name4)
75
76print 'add G1 + G2'
77G5 = G1+G2
78
79#-------------------------------------------------------------------------------
80# Combine, clip and export dataset
81#-------------------------------------------------------------------------------
82
83print 'clipping outside extents'
84#G1_clip = G1.clip(project.extent_elev_dir4, verbose=True)
85G3_clip1 = G3.clip_outside(project.extent_elev_dir1, verbose=True) 
86G3_clip2 = G3_clip1.clip_outside(project.extent_elev_dir2, verbose=True)
87G4_clip = G4.clip_outside(project.extent_elev_dir3, verbose=True)
88
89#print 'hello', dir(G3_clip1), dir(G3_clip1), dir(G4_clip)
90print 'Add all geospatial objects'
91G = G1 + G2 + G3_clip2 + G4_clip
92
93print 'Clip combined geospatial object by bounding polygon'
94G_clipped = G.clip(project.poly_all)
95
96print 'Oxport combined DEM file'
97G_clipped.export_points_file(project.combined_dir_name + '.pts')
98G_clipped.export_points_file(project.combined_dir_name + '.txt') #Use for comparision in ARC
99
Note: See TracBrowser for help on using the repository browser.