source: anuga_work/production/dampier_2006/project.py @ 3997

Last change on this file since 3997 was 3997, checked in by nick, 18 years ago

updates to dampier

File size: 12.5 KB
Line 
1"""Common filenames and locations for topographic data, meshes and outputs.
2"""
3
4import sys
5from os import sep, environ, getenv, getcwd
6from os.path import expanduser, basename
7
8from anuga.coordinate_transforms.redfearn import\
9     degminsec2decimal_degrees,\
10     convert_from_latlon_to_utm
11
12from time import localtime, strftime, gmtime, ctime
13from anuga.geospatial_data.geospatial_data import *
14from anuga.utilities.polygon import read_polygon, plot_polygons, polygon_area
15from anuga.utilities.system_tools import get_user_name
16
17# file and system info
18#---------------------------------
19codename = 'project.py'
20
21home = getenv('INUNDATIONHOME') #Sandpit's parent dir   
22user = get_user_name()
23
24# INUNDATIONHOME is the inundation directory, not the data directory.
25home += sep +'data'
26#----------------------------------
27# Location and naming of scenario data
28#----------------------------------
29state = 'western_australia'
30scenario_name = 'dampier_tsunami'
31scenario_datas_name = 'dampier_tsunami_scenario_2006'  #name of the directory where the data is stored
32#scenario_datas_name = 'karratha_tsunami_scenario_2005' # Tmp location
33
34#mesh_name = 'elevation50m'
35boundaries_name = 'dampier'
36boundaries_source = 'mag_9_corrected'
37#boundaries_source = 'test'
38
39tide = 2.4
40#tide = 0.0
41
42# topography file names
43onshore_name = 'dli_no_islands'
44coast_name = 'DTED_05_Contour'
45islands_name = 'dted_islands'
46offshore_name = 'XY100003902'
47offshore_name1 = 'XY100003903'
48offshore_name2 = 'XY100003951'
49offshore_name3 = 'XY100006321'
50offshore_name4 = 'XY100011756'
51offshore_name5 = 'XY100014243'
52offshore_name6 = 'XY100014244'
53offshore_name7 = 'XY100021081'
54offshore_name8 = 'XY100021082'
55offshore_name9 = 'XY100021083'
56offshore_name10 = 'XY100021085'
57offshore_name11 = 'XY100021086'
58offshore_name12 = 'XY100026309'
59offshore_name13 = 'XY100026338'
60offshore_name14 = 'XYDM83'
61
62offshore_old = 'elevation50m'
63
64combined_name ='dampier_combined_elevation'
65combined_final_name ='dampier_combined_elevation_final'
66
67gauge_name = 'dampier_gauges_up2.csv'
68
69#Derive subdirectories and filenames
70
71meshes_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'meshes'+sep
72topographies_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'topographies'+sep
73gauges_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'gauges'+sep
74polygons_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'polygons'+sep
75boundaries_in_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'boundaries'+sep+'urs'+sep+boundaries_source+sep
76#outputdir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'output'+sep
77tide_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'tide_data'+sep
78
79time = strftime('%Y%m%d_%H%M%S',localtime()) #gets time for new dir
80gtime = strftime('%Y%m%d_%H%M%S',gmtime()) #gets time for new dir
81#cctime = strftime('%Y%m%d_%H%M%S',ctime()) #gets time for new dir
82build_time = time+'_build'
83run_time = time+'_run'
84
85print 'gtime: ', gtime
86output_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'outputs'+sep
87output_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'outputs'+sep
88output_build_time_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'outputs'+sep+build_time+sep
89output_run_time_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'outputs'+sep+run_time+sep
90topographies_time_dir = topographies_dir+build_time+sep
91boundaries_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'boundaries'+sep
92boundaries_time_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'boundaries'+sep+build_time+sep
93meshes_time_dir = meshes_dir+build_time+sep
94
95#ideas
96#boundaries_time_dir = boundaries_in_dir+'urs'+sep+boundaries_source+sep
97
98gauge_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'gauges'+sep
99gauge_filename = gauge_dir + 'dampier_gauges_up2.csv'
100
101gauges_dir_name = gauges_dir + gauge_name
102
103onshore_dir_name = topographies_dir + onshore_name
104coast_dir_name = topographies_dir + coast_name
105islands_dir_name = topographies_dir + islands_name
106offshore_dir_name = topographies_dir + offshore_name
107offshore_dir_name1 = topographies_dir + offshore_name1
108offshore_dir_name2 = topographies_dir + offshore_name2
109offshore_dir_name3 = topographies_dir + offshore_name3
110offshore_dir_name4 = topographies_dir + offshore_name4
111offshore_dir_name5 = topographies_dir + offshore_name5
112offshore_dir_name6 = topographies_dir + offshore_name6
113offshore_dir_name7 = topographies_dir + offshore_name7
114offshore_dir_name8 = topographies_dir + offshore_name8
115offshore_dir_name9 = topographies_dir + offshore_name9
116offshore_dir_name10 = topographies_dir + offshore_name10
117offshore_dir_name11 = topographies_dir + offshore_name11
118offshore_dir_name12 = topographies_dir + offshore_name12
119offshore_dir_name13 = topographies_dir + offshore_name13
120offshore_dir_name14 = topographies_dir + offshore_name14
121
122offshore_dir_name_old = topographies_dir + offshore_old
123
124
125
126#output dir
127combined_dir_name = topographies_dir + combined_name
128combined_time_dir_name = topographies_time_dir + combined_name
129combined_time_dir_final_name = topographies_time_dir + combined_final_name
130
131
132meshes_dir_name = meshes_dir + scenario_name
133meshes_time_dir_name = meshes_time_dir + scenario_name
134#output_build_time_dir_name = output_build_time_dir + scenario_name  #Used by post processing
135output_run_time_dir_name = output_run_time_dir + scenario_name  #Used by post processing
136boundaries_in_dir_name = boundaries_in_dir + boundaries_name
137boundaries_time_dir_name = boundaries_time_dir + boundaries_name  #Used by post processing
138boundaries_dir_name = boundaries_dir + boundaries_name
139
140
141
142# Regions
143
144refzone = 50 
145south = degminsec2decimal_degrees(-20,55,0)
146north = degminsec2decimal_degrees(-20,15,0)
147#north = degminsec2decimal_degrees(-19,15,0)
148west = degminsec2decimal_degrees(116,17,0)
149east = degminsec2decimal_degrees(117,10,0)
150
151#only used to clip boundary condition
152#north_boundary = north + 0.02
153#south_boundary = south - 0.02
154#west_boundary = west - 0.02
155#east_boundary = east + 0.02
156
157south_boundary = degminsec2decimal_degrees(-21,0,0)
158#north_boundary = degminsec2decimal_degrees(-19,00,0)
159north_boundary = degminsec2decimal_degrees(-20,10,0)
160#west_boundary = degminsec2decimal_degrees(116,0,0)
161#east_boundary = degminsec2decimal_degrees(118,00,0)
162west_boundary = degminsec2decimal_degrees(116,10,0)
163east_boundary = degminsec2decimal_degrees(117,15,0)
164
165
166p0 = [south, degminsec2decimal_degrees(116,32,0)]
167p1 = [south, west]
168p2 = [degminsec2decimal_degrees(-20,23,0), west]
169p3 = [north, degminsec2decimal_degrees(116,45,0)]
170p4 = [north, degminsec2decimal_degrees(117,0,0)]
171p5 = [p2[0], degminsec2decimal_degrees(117,8,0)]
172p6 = [degminsec2decimal_degrees(-20,30,0), east]
173p7 = [degminsec2decimal_degrees(-20,38,0), east]
174p8 = [south, east]
175
176bounding_polygon, zone =\
177                  convert_from_latlon_to_utm([p0, p1, p2, p3, p4, p5, p6, p7, p8])
178#bounding_polygon, zone =\
179#                  convert_from_latlon_to_utm([p1, p2, p3, p4, p5, p6, p7])
180refzone = zone
181from anuga.utilities.polygon import read_polygon, plot_polygons, polygon_area, is_inside_polygon
182print 'poly area', polygon_area(bounding_polygon)/1000000.0
183
184#Interior regions
185
186# CIPMA point of interest
187cipma_latitude = -20.588456
188cipma_longitude = 116.771527
189
190
191k0 = [cipma_latitude-0.02, cipma_longitude-0.02]
192k1 = [cipma_latitude-0.02, cipma_longitude+0.02]
193k2 = [cipma_latitude+0.02, cipma_longitude+0.02]
194k3 = [cipma_latitude+0.02, cipma_longitude-0.02]
195
196cipma_polygon, zone = convert_from_latlon_to_utm([k0, k1, k2, k3])
197assert zone == refzone
198
199e_min_area = 474000
200e_max_area = 480000
201n_min_area = 7719000
202n_max_area = 7725000
203
204poly_facility = read_polygon(polygons_dir+'facility.csv')
205
206poly_pipeline = read_polygon(polygons_dir+'pipeline2.csv')
207
208poly_interior = read_polygon(polygons_dir+'interior.csv')
209
210poly_coast = read_polygon(polygons_dir+'coast_final.csv')
211
212clip_poly_e = read_polygon(polygons_dir+'gap_e.csv')
213
214clip_poly_nw = read_polygon(polygons_dir+'gap_nw.csv')
215
216clip_poly_mid_w = read_polygon(polygons_dir+'gap_mid_w.cvs')
217
218clip_poly_mid_e = read_polygon(polygons_dir+'gap_mid_e.cvs')
219
220#Interior regions
221karratha_south = degminsec2decimal_degrees(-20,44,0)
222karratha_north = degminsec2decimal_degrees(-20,42,0)
223karratha_west = degminsec2decimal_degrees(116,48,0)
224karratha_east = degminsec2decimal_degrees(116,53,30)
225
226k0 = [karratha_south, karratha_west]
227k1 = [karratha_south, karratha_east]
228k2 = [karratha_north, karratha_east]
229k3 = [karratha_north, karratha_west]   
230
231karratha_polygon, zone = convert_from_latlon_to_utm([k0, k1, k2, k3])
232assert zone == refzone
233
234
235#Interior regions
236dampier_south = degminsec2decimal_degrees(-20,40,0)
237dampier_north = degminsec2decimal_degrees(-20,38,10)
238dampier_west = degminsec2decimal_degrees(116,43,0)
239dampier_east = degminsec2decimal_degrees(116,45,0)
240
241d0 = [dampier_south, dampier_west]
242d1 = [dampier_south, dampier_east]
243d2 = [dampier_north, dampier_east]
244d3 = [dampier_north, dampier_west]   
245
246dampier_polygon, zone = convert_from_latlon_to_utm([d0, d1, d2, d3])
247assert zone == refzone
248
249
250#Interior regions
251refinery_south = degminsec2decimal_degrees(-20,37,50)
252refinery_north = degminsec2decimal_degrees(-20,36,0)
253refinery_west = degminsec2decimal_degrees(116,44,0)
254refinery_east = degminsec2decimal_degrees(116,46,10)
255
256d0 = [refinery_south, refinery_west]
257d1 = [refinery_south, refinery_east]
258d2 = [refinery_north, refinery_east]
259d3 = [refinery_north, refinery_west]   
260
261refinery_polygon, zone = convert_from_latlon_to_utm([d0, d1, d2, d3])
262assert zone == refzone
263
264
265#Interior region around 468899, 7715177:
266#lat (-20, 39, 44.93753), lon (116, 42, 5.09106)
267
268point_south = degminsec2decimal_degrees(-20,39,46)
269point_north = degminsec2decimal_degrees(-20,39,42)
270point_west = degminsec2decimal_degrees(116,42,0)
271point_east = degminsec2decimal_degrees(116,42,10)
272
273d0 = [point_south, point_west]
274d1 = [point_south, point_east]
275d2 = [point_north, point_east]
276d3 = [point_north, point_west]   
277
278point_polygon, zone = convert_from_latlon_to_utm([d0, d1, d2, d3])
279assert zone == refzone
280
281
282#Neils areas around interesting points
283neil1_point1 = [degminsec2decimal_degrees(-20,35,34),
284                degminsec2decimal_degrees(116,45,18)]
285neil1_point2 = [degminsec2decimal_degrees(-20,36,15),
286                degminsec2decimal_degrees(116,46,18)]
287neil1_point3 = [degminsec2decimal_degrees(-20,35,9),
288                degminsec2decimal_degrees(116,47,17)]
289neil1_point4 = [degminsec2decimal_degrees(-20,34,26),
290                degminsec2decimal_degrees(116,46,17)]
291
292neil1_polygon, zone = convert_from_latlon_to_utm([neil1_point1,
293                                                  neil1_point2,
294                                                  neil1_point3,
295                                                  neil1_point4])
296assert zone == refzone
297
298
299
300
301neil2_point1 = [degminsec2decimal_degrees(-20,39,36),
302                degminsec2decimal_degrees(116,41,33)]
303neil2_point2 = [degminsec2decimal_degrees(-20,40,10),
304                degminsec2decimal_degrees(116,42,13)]
305neil2_point3 = [degminsec2decimal_degrees(-20,38,39),
306                degminsec2decimal_degrees(116,43,49)]
307neil2_point4 = [degminsec2decimal_degrees(-20,38,5),
308                degminsec2decimal_degrees(116,43,9)]
309
310neil2_polygon, zone = convert_from_latlon_to_utm([neil2_point1,
311                                                  neil2_point2,
312                                                  neil2_point3,
313                                                  neil2_point4])
314assert zone == refzone
315
316
317
318
319
320#Withnell bay
321wb_point1 = [degminsec2decimal_degrees(-20,35,34),
322                degminsec2decimal_degrees(116,45,18)]
323wb_point2 = [degminsec2decimal_degrees(-20,36,15),
324                degminsec2decimal_degrees(116,46,18)]
325wb_point3 = [degminsec2decimal_degrees(-20,35,9),
326                degminsec2decimal_degrees(116,47,17)]
327wb_point4 = [degminsec2decimal_degrees(-20,34,26),
328                degminsec2decimal_degrees(116,46,17)]
329
330wb_polygon, zone = convert_from_latlon_to_utm([wb_point1, wb_point2,
331                                               wb_point3, wb_point4])
332assert zone == refzone
333
334
335
336
337
338#Larger Withnell bay
339lwb_point1 = [degminsec2decimal_degrees(-20,35,59),
340                degminsec2decimal_degrees(116,42,00)]
341lwb_point2 = [degminsec2decimal_degrees(-20,36,50),
342                degminsec2decimal_degrees(116,46,50)]
343lwb_point3 = [degminsec2decimal_degrees(-20,34,00),
344                degminsec2decimal_degrees(116,47,39)]
345lwb_point4 = [degminsec2decimal_degrees(-20,33,00),
346                degminsec2decimal_degrees(116,42,50)]
347
348lwb_polygon, zone = convert_from_latlon_to_utm([lwb_point1, lwb_point2,
349                                                lwb_point3, lwb_point4])
350                                                     
351assert zone == refzone
352
353
354
355
Note: See TracBrowser for help on using the repository browser.