source: anuga_work/production/shark_bay_2007/project.py @ 4624

Last change on this file since 4624 was 4624, checked in by sexton, 17 years ago

export scripts for Shark Bay

File size: 7.4 KB
Line 
1"""Common filenames and locations for topographic data, meshes and outputs.
2Also includes origin for slump scenario.
3
4All filenames are given without extension
5"""
6
7from os import sep, environ, getenv, getcwd, umask
8from os.path import expanduser, basename, join
9from anuga.utilities.polygon import read_polygon, plot_polygons, polygon_area, is_inside_polygon, number_mesh_triangles
10import sys
11from anuga.coordinate_transforms.redfearn import degminsec2decimal_degrees
12from time import localtime, strftime, gmtime
13from anuga.utilities.system_tools import get_user_name, get_host_name
14
15codename = 'project.py' # FIXME can be obtained automatically as __file__
16
17home = join(getenv('INUNDATIONHOME'), 'data') # Location of Data   
18user = get_user_name()
19host = get_host_name()
20#needed when running using mpirun, mpirun doesn't inherit umask from .bashrc
21umask(002)
22
23#Making assumptions about the location of scenario data
24state = 'western_australia'
25scenario_name = 'shark_bay'
26scenario = 'shark_bay_tsunami_scenario'
27
28#time stuff
29time = strftime('%Y%m%d_%H%M%S',gmtime()) #gets time for new dir
30build_time = time+'_build'
31run_time = time+'_run'
32
33#tide = 0.0 # MSL
34tide = 0.85 # HAT average between Denham and Canarvon
35
36#Maybe will try to make project a class to allow these parameters to be passed in.
37alpha = 0.1
38friction=0.01
39finaltime=40000
40starttime=6000 # Action starts at 9000
41setup='final'
42#setup='trial'
43source='shark_bay'
44
45#boundary_event = 'july2006'
46boundary_event = '10000'
47
48if setup =='trial':
49    print'trial'
50    res_factor=10
51    time_thinning=48
52    yieldstep=240
53if setup =='basic': 
54    print'basic'
55    res_factor=4
56    time_thinning=12
57    yieldstep=120
58if setup =='final': 
59    print'final'
60    res_factor=1.0
61    time_thinning=4
62    yieldstep=60
63
64
65
66dir_comment='_'+setup+'_'+str(tide)+'_'+str(source)+'_'+\
67             str(user)+'_'+boundary_event+'_'+'tide'+str(tide)+'_coastpolys'
68
69#dir_comment='_'+setup+'_'+str(tide)+'_'+str(source)+'_'+\
70#             str(user)+'_'+boundary_event
71
72# elevation data filenames
73ascii_grid_filenames = ['10m_dem_without_survey', '50m_dem_without_10m_dem']
74point_filenames = ['field_survey_north', 'field_survey_south',
75                   'clipped_bathymetry_final', 'coast_points_final']
76
77# final topo name
78combined_name ='shark_bay_combined_elevation'
79combined_small_name = 'shark_bay_combined_elevation_small'
80
81
82anuga_dir = join(home, state, scenario, 'anuga')
83
84topographies_in_dir = join(home, state, scenario,
85                           'elevation_final', 'test_area')
86
87topographies_dir = join(anuga_dir, 'topographies') # Output dir for ANUGA
88
89# Convert topo file locations into absolute pathnames
90for filename_list in [ascii_grid_filenames, point_filenames]:
91    for i, filename in enumerate(filename_list):
92        filename_list[i] = join(topographies_in_dir, filename)
93
94#final topo files
95combined_dir_name = join(topographies_dir, combined_name)
96combined_small_dir_name = join(topographies_dir, combined_small_name)
97
98
99mesh_dir = join(anuga_dir, 'meshes')
100mesh_name = join(mesh_dir, scenario_name)
101
102
103polygons_dir = join(anuga_dir, 'polygons') # Created with ArcGIS (csv files)
104tide_dir = join(anuga_dir, 'tide_data')
105
106
107# locations for boundary conditions
108if source=='shark_bay':
109    if boundary_event == '10000':
110        boundary_file_name = 'shark_bay_3867_18052007'
111        boundary_dir = join(anuga_dir, 'boundaries', 'shark_bay', '1_10000')
112    elif boundary_event == 'july2006':
113        boundary_file_name = 'july2006'   
114        boundary_dir = join(anuga_dir, 'boundaries',
115                            'shark_bay', 'july_2006_event')
116    else:
117        raise Exception, 'Unknown boundary event specified'
118
119boundary_name = join(boundary_dir, boundary_file_name)
120
121#output locations
122output_dir = join(anuga_dir, 'outputs')+sep
123output_build_time_dir = output_dir+build_time+sep
124output_run_time_dir = output_dir +run_time+dir_comment+sep
125output_run_time_dir_name = output_run_time_dir + scenario_name  #Used by post processing
126
127#gauges
128gauge_name = 'gaugesv2.csv' 
129#gauge_name_test = 'gauge_checking_test.csv'
130gauges_dir = anuga_dir+sep+'gauges'+sep
131gauges_dir_name = gauges_dir + gauge_name
132#gauges_dir_name_test = gauges_dir + gauge_name_test
133
134#tide_dir = anuga_dir+'tide_data'+sep
135
136#buildings_filename = gauges_dir + 'pt_hedland_res.csv'
137#buildings_filename_out = 'pt_hedland_res_modified.csv'
138#buildings_filename_damage_out = 'pt_hedland_res_modified_damage.csv'
139#tidal_filename = tide_dir + 'pt_hedland_tide.txt'
140#community_filename = gauges_dir + 'CHINS_v2.csv'
141#community_scenario = gauges_dir + 'community_pt_hedland.csv'
142
143# clipping region to make DEM (pts file) from onshore data
144#eastingmin = 594000
145#eastingmax = 715000
146#northingmin = 7720000
147#northingmax = 7880000
148
149# for ferret2sww
150south = degminsec2decimal_degrees(-20,30,0)
151north = degminsec2decimal_degrees(-17,10,0)
152west = degminsec2decimal_degrees(117,00,0)
153east = degminsec2decimal_degrees(120,00,0)
154
155# region to export (used from export_results.py)
156e_min_area = 713500   # Eastings min
157e_max_area = 726350
158n_min_area = 7100890
159n_max_area = 7111250
160
161#export_region = [[e_min_area,n_min_area],[e_min_area,n_max_area],
162#                 [e_max_area,n_max_area],[e_max_area,n_min_area]]
163                 
164#FIXME (Ole): It is bad to read in polygons in project.py. If a file does not exist
165#project.py cannot be imported
166
167# Bounding polygon
168bounding_polygon = read_polygon(join(polygons_dir, 'boundary_extent_small.csv'))
169res_bounding_polygon = 500000*res_factor
170
171#Interior regions
172interior_regions_and_resolutions = {
173    'study_area': 300000,
174    'bernier_dorre_islands': 250000,
175    'map_area': 10000,
176    'bay_area': 4000,
177    'south_bank_buffer': 2000,
178    'south_coast_polygon_1': 1000,
179    'south_coast_polygon_2': 1000,
180    'north_bay_coastline': 1000,
181    'tsunami_approach_area': 50,
182    'inundated_area': 10,
183    'exclusion_1': 10000000,
184    'exclusion_2': 10000000,
185    'exclusion_3': 10000000}
186
187
188#poly_exclusion_area1 = read_polygon(join(polygons_dir, 'exclusion_area1.csv'))
189#res_exclusion_area1 = 2000000*res_factor
190#
191#interior_regions = [[poly_inundated_area, res_inundated_area],
192#                    [poly_tsunami_approach_area, res_tsunami_approach_area],
193#                    [poly_bay_area, res_bay_area],                   
194#                    [poly_map_area, res_map_area]]
195                    #[poly_model_area, res_model_area],
196                    #[poly_exclusion_area1, res_exclusion_area1]]                                       
197
198# FIXME (Ole): What if tags are wrong?
199
200#Large bounding poly
201#boundary_tags = {'back': [1, 2, 3, 4, 5],   
202#                 'side': [0, 6],
203#                 'ocean': [7, 8, 9, 10, 11]}
204
205
206#Small bounding poly
207#(counter clockwise starting in upper right corner)
208boundary_tags = {'tide': [0, 1, 2],
209                 'ocean': [3, 4, 5, 6]}
210
211
212
213# Read filenames etc (FIXME: move this to new file: read_structures)
214
215interior_regions = [] # Consider using dictionary for mesh interface.
216for area_name in interior_regions_and_resolutions.keys():
217    file_name = join(polygons_dir, area_name + '.csv')
218   
219    print 'Reading polygon from', file_name
220    polygon = read_polygon(file_name)
221    base_resolution = interior_regions_and_resolutions[area_name]
222   
223    interior_regions.append([polygon, base_resolution*res_factor])
224
225
226
227trigs_min = number_mesh_triangles(interior_regions,
228                                  bounding_polygon,
229                                  res_bounding_polygon)
230
231
232onshore_polygon = read_polygon(join(topographies_in_dir,
233                                    'initial_condition.txt'))
234
Note: See TracBrowser for help on using the repository browser.