1 | """ |
---|
2 | Script for building the elevation data to run a tsunami inundation scenario |
---|
3 | for patong, Thailand. |
---|
4 | |
---|
5 | Input: elevation data from project.py |
---|
6 | Output: pts file stored in project.topographies_dir |
---|
7 | The run_patong.py is reliant on the output of this script. |
---|
8 | |
---|
9 | """ |
---|
10 | |
---|
11 | #------------------------------------------------------------------------------ |
---|
12 | # Import necessary modules |
---|
13 | #------------------------------------------------------------------------------ |
---|
14 | |
---|
15 | # Standard modules |
---|
16 | from os import sep |
---|
17 | from os.path import dirname, basename |
---|
18 | from os import mkdir, access, F_OK |
---|
19 | from shutil import copy |
---|
20 | import time |
---|
21 | import sys |
---|
22 | |
---|
23 | # Related major packages |
---|
24 | from anuga.shallow_water import Domain |
---|
25 | from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf, dem2pts |
---|
26 | from anuga.geospatial_data.geospatial_data import * |
---|
27 | from anuga.shallow_water.data_manager import start_screen_catcher, copy_code_files |
---|
28 | |
---|
29 | # Application specific imports |
---|
30 | import 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 | |
---|
37 | copy_code_files(project.output_build_time_dir,__file__, |
---|
38 | dirname(project.__file__)+sep+ project.__name__+'.py' ) |
---|
39 | |
---|
40 | start_screen_catcher(project.output_build_time_dir) |
---|
41 | |
---|
42 | print '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 | #------------------------------------------------------------------------------- |
---|
51 | print "project.poly_all", project.poly_all |
---|
52 | print "project.combined_dir_name", project.combined_dir_name |
---|
53 | |
---|
54 | # input elevation directory filenames |
---|
55 | elevation_in_dir_name1 = project.elevation_in_dir_name1 |
---|
56 | elevation_in_dir_name2 = project.elevation_in_dir_name2 |
---|
57 | elevation_in_dir_name3 = project.elevation_in_dir_name3 |
---|
58 | elevation_in_dir_name4 = project.elevation_in_dir_name4 |
---|
59 | |
---|
60 | # output elevation directory filenames |
---|
61 | elevation_dir_name1 = project.elevation_dir_name1 |
---|
62 | elevation_dir_name2 = project.elevation_dir_name2 |
---|
63 | elevation_dir_name3 = project.elevation_dir_name3 |
---|
64 | elevation_dir_name4 = project.elevation_dir_name4 |
---|
65 | |
---|
66 | # create elevation pts files |
---|
67 | print'create Geospatial data1 objects from topographies' |
---|
68 | G1 = Geospatial_data(file_name = elevation_in_dir_name1) |
---|
69 | print'create Geospatial data2 objects from topographies' |
---|
70 | G2 = Geospatial_data(file_name = elevation_in_dir_name2) |
---|
71 | print'create Geospatial data3 objects from topographies' |
---|
72 | G3 = Geospatial_data(file_name = elevation_in_dir_name3) |
---|
73 | print'create Geospatial data4 objects from topographies' |
---|
74 | G4 = Geospatial_data(file_name = elevation_in_dir_name4) |
---|
75 | |
---|
76 | print 'add G1 + G2' |
---|
77 | G5 = G1+G2 |
---|
78 | |
---|
79 | #------------------------------------------------------------------------------- |
---|
80 | # Combine, clip and export dataset |
---|
81 | #------------------------------------------------------------------------------- |
---|
82 | |
---|
83 | print 'clipping outside extents' |
---|
84 | #G1_clip = G1.clip(project.extent_elev_dir4, verbose=True) |
---|
85 | G3_clip1 = G3.clip_outside(project.extent_elev_dir1, verbose=True) |
---|
86 | G3_clip2 = G3_clip1.clip_outside(project.extent_elev_dir2, verbose=True) |
---|
87 | G4_clip = G4.clip_outside(project.extent_elev_dir3, verbose=True) |
---|
88 | |
---|
89 | #print 'hello', dir(G3_clip1), dir(G3_clip1), dir(G4_clip) |
---|
90 | print 'Add all geospatial objects' |
---|
91 | G = G1 + G2 + G3_clip2 + G4_clip |
---|
92 | |
---|
93 | print 'Clip combined geospatial object by bounding polygon' |
---|
94 | G_clipped = G.clip(project.poly_all) |
---|
95 | |
---|
96 | print 'Oxport combined DEM file' |
---|
97 | G_clipped.export_points_file(project.combined_dir_name + '.pts') |
---|
98 | G_clipped.export_points_file(project.combined_dir_name + '.txt') #Use for comparision in ARC |
---|
99 | |
---|