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

Last change on this file since 4311 was 4311, checked in by nick, 17 years ago

updates to dampier
File size: 15.7 KB
Line 
1"""Common filenames and locations for topographic data, meshes and outputs.
2"""
3
4from os import sep, environ, getenv, getcwd
5from os.path import expanduser
6import sys
7from time import localtime, strftime, gmtime
8from anuga.utilities.polygon import read_polygon, plot_polygons, polygon_area, is_inside_polygon, number_mesh_triangles
9from anuga.coordinate_transforms.redfearn import degminsec2decimal_degrees, convert_from_latlon_to_utm
10from anuga.utilities.system_tools import get_user_name
11
12# file and system info
13#---------------------------------
14codename = 'project.py'
15
16home = getenv('INUNDATIONHOME') #Sandpit's parent dir   
17user = get_user_name()
18
19# INUNDATIONHOME is the inundation directory, not the data directory.
20home += sep +'data'
21
22#time stuff
23time = strftime('%Y%m%d_%H%M%S',localtime()) #gets time for new dir
24gtime = strftime('%Y%m%d_%H%M%S',gmtime()) #gets time for new dir
25build_time = time+'_build'
26run_time = time+'_run'
27temp_time = time+'_temp'
28print 'gtime: ', gtime
29
30tide = 2.4
31
32#Making assumptions about the location of scenario data
33state = 'western_australia'
34scenario_name = 'dampier'
35scenario = 'dampier_tsunami_scenario_2006'
36
37# onshore data provided by WA DLI
38onshore_name = 'DLI_DTED_raster_clipped' # original
39
40# AHO + DPI data + colin French coastline
41coast_name = 'coastline_edited_w_DEM'
42offshore_name = 'clipped_bathy'
43offshore1_name = 'elev_501'
44offshore2_name = 'inferrec_e'
45
46#final topo name
47combined_name ='dampier_combined_elevation'
48combined_smaller_name = 'dampier_combined_elevation_small'
49combined_smallest_name = 'dampier_combined_elevation_smallest'
50
51topographies_in_dir = home+sep+state+sep+scenario+sep+'elevation_final'+sep+'070112'+sep+'points'+sep
52topographies_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'topographies'+sep
53topographies_time_dir = topographies_dir+build_time+sep
54
55temp_dir = home+sep+state+sep+scenario+sep+'anuga'+sep
56temp_dir_name = temp_dir + temp_time+sep
57
58# input topo file location
59onshore_in_dir_name = topographies_in_dir + onshore_name
60coast_in_dir_name = topographies_in_dir + coast_name
61offshore_in_dir_name = topographies_in_dir + offshore_name
62offshore1_in_dir_name = topographies_in_dir + offshore1_name
63offshore2_in_dir_name = topographies_in_dir + offshore2_name
64
65onshore_dir_name = topographies_dir + onshore_name
66coast_dir_name = topographies_dir + coast_name
67offshore_dir_name = topographies_dir + offshore_name
68offshore1_dir_name = topographies_dir + offshore1_name
69offshore2_dir_name = topographies_dir + offshore2_name
70
71#final topo files
72combined_dir_name = topographies_dir + combined_name
73combined_time_dir_name = topographies_time_dir + combined_name
74combined_smaller_dir_name = topographies_dir + combined_smaller_name
75combined_smallest_dir_name = topographies_dir + combined_smallest_name
76#combined_time_dir_final_name = topographies_time_dir + combined_final_name
77
78meshes_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'meshes'+sep
79meshes_dir_name = meshes_dir + scenario_name
80
81polygons_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'polygons'+sep+'2007polys'+sep
82tide_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'tide_data'+sep
83
84boundaries_source = ''
85boundaries_name = 'o'
86boundaries_name1 = 'o_new1'
87boundaries_name2 = 'o_test'
88
89#boundaries locations
90#boundaries_in_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'boundaries'+sep
91boundaries_in_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'boundaries'+sep+'urs'+sep+'1_10000'+sep
92boundaries_in_dir_name = boundaries_in_dir + boundaries_name
93boundaries_in_dir2 = home+sep+state+sep+scenario+sep+'anuga'+sep+'boundaries'+sep+'urs'+sep+'ungridded'+sep
94boundaries_in_dir_name2 = boundaries_in_dir2 + boundaries_name
95boundaries_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'boundaries'+sep
96boundaries_dir_name = boundaries_dir + boundaries_name
97boundaries_dir_name1 = boundaries_dir + boundaries_name1
98boundaries_dir_name2 = boundaries_dir + boundaries_name2
99boundaries_dir_name3 = boundaries_dir + boundaries_name+'_test_8500_12000'
100boundaries_dir_name4 = boundaries_dir + boundaries_name+'_8500_12000_no_zone'
101boundaries_dir_name5 = boundaries_dir + '20061026_173920_build'+sep+'dampier'
102boundaries_dir_name6 = boundaries_dir + 'SU-AU'
103
104#boundaries_time_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'boundaries'+sep+build_time+sep
105#boundaries_time_dir_name = boundaries_time_dir + boundaries_name  #Used by post processing
106
107#output locations
108output_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'outputs'+sep
109output_build_time_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'outputs'+sep+build_time+sep
110output_run_time_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'outputs'+sep+run_time+sep
111output_run_time_dir_name = output_run_time_dir + scenario_name  #Used by post processing
112
113#gauges
114gauge_name = 'dampier_gauges_up2.csv' #'dampier.csv'
115gauge_name_simple = 'dampier_gauges_simple.csv' 
116gauge_name_test = 'dampier_gauges_up_test.csv' #'dampier.csv'
117gauges_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'gauges'+sep
118gauges_dir_name = gauges_dir + gauge_name
119gauges_dir_name_simple = gauges_dir + gauge_name_simple
120gauges_dir_name_test = gauges_dir + gauge_name_test
121
122
123#buildings_filename = gauges_dir + 'dampier_res_nexis.csv'
124buildings_filename_damage_out = 'dampier_res_nexis_modified.csv'
125###############################
126# Domain definitions
127###############################
128
129refzone = 50 
130south_boundary = degminsec2decimal_degrees(-20,58,0)
131#north_boundary = degminsec2decimal_degrees(-20,13,0)
132north_boundary = degminsec2decimal_degrees(-19,30,0)
133west_boundary = degminsec2decimal_degrees(116,15,0)
134east_boundary = degminsec2decimal_degrees(117,11,0)
135##
136##p0 = [south, degminsec2decimal_degrees(116,32,0)]
137##p1 = [south, west]
138##p2 = [degminsec2decimal_degrees(-20,23,0), west]
139##p3 = [north, degminsec2decimal_degrees(116,45,0)]
140##p4 = [north, degminsec2decimal_degrees(117,0,0)]
141##p5 = [p2[0], degminsec2decimal_degrees(117,8,0)]
142##p6 = [degminsec2decimal_degrees(-20,30,0), east]
143##p7 = [degminsec2decimal_degrees(-20,38,0), east]
144##p8 = [south, east]
145##
146##poly_all, zone = convert_from_latlon_to_utm([p0, p1, p2, p3, p4, p5, p6, p7, p8])
147##refzone = zone
148poly_all = read_polygon(polygons_dir+'extent.csv')
149print 'Area of bounding polygon', polygon_area(poly_all)/1000000.0
150
151res_poly_all = 400000
152#res_poly_all = 500000
153
154###############################
155# Interior region definitions
156###############################
157
158poly_region = read_polygon(polygons_dir+'region.csv')
159res_region = 50000
160
161poly_dampier = read_polygon(polygons_dir+'dampier_town.csv')
162res_dampier = 1000
163#res_dampier = 500
164#res_dampier = 5000
165
166poly_karratha = read_polygon(polygons_dir+'karrathav2.csv')
167res_karratha = 15000
168
169poly_karratha_town = read_polygon(polygons_dir+'karratha_townv2.csv')
170res_karratha_town = 1000
171#res_karratha_town = 500
172#res_karratha_town = 5000
173
174poly_facility = read_polygon(polygons_dir+'facility.csv')
175res_facility = 1000
176
177poly_delambre = read_polygon(polygons_dir+'delambre.csv')
178res_delambre = 5000
179
180poly_coast = read_polygon(polygons_dir+'coastpoly.csv')
181#res_coast = 1000
182res_coast = 10000
183
184poly_NWislands = read_polygon(polygons_dir+'nw_islands_area.csv')
185res_NWislands = 50000
186
187poly_island0 = read_polygon(polygons_dir+'island0.csv')
188res_island0 = res_poly_all
189
190poly_island1 = read_polygon(polygons_dir+'island1.csv')
191res_island0 = res_poly_all
192
193poly_island2 = read_polygon(polygons_dir+'island2.csv')
194res_island0 = res_poly_all
195
196poly_island3 = read_polygon(polygons_dir+'island3.csv')
197res_island0 = res_poly_all
198
199res_islands = 5000
200#res_islands = 15000
201
202poly_ref_nw4 = read_polygon(polygons_dir+'ref_nw4.csv')
203res_ref_nw4 = res_islands
204
205poly_island4 = read_polygon(polygons_dir+'island4.csv')
206res_island0 = res_poly_all
207
208poly_ref_nw5 = read_polygon(polygons_dir+'ref_nw5.csv')
209res_ref_nw5 = res_islands
210
211poly_island5 = read_polygon(polygons_dir+'island5.csv')
212res_island0 = res_poly_all
213
214poly_ref_nw6 = read_polygon(polygons_dir+'ref_nw6.csv')
215res_ref_nw6 = res_islands
216
217poly_island6 = read_polygon(polygons_dir+'island6.csv')
218res_island0 = res_poly_all
219
220poly_ref_nw7 = read_polygon(polygons_dir+'ref_nw7.csv')
221res_ref_nw7 = res_islands
222
223poly_island7 = read_polygon(polygons_dir+'island7.csv')
224res_island0 = res_poly_all
225
226poly_ref_nw8 = read_polygon(polygons_dir+'ref_nw8.csv')
227res_ref_nw8 = res_islands
228
229poly_island8 = read_polygon(polygons_dir+'island8.csv')
230res_island0 = res_poly_all
231
232
233##plot_polygons([poly_dampier,poly_karratha,poly_karratha_town,poly_delambre,
234##                poly_coast,poly_NWislands,poly_island0,poly_island1,poly_island2,
235##                poly_island3,poly_island4,poly_island5,poly_island6,
236##                poly_island7,poly_island8,poly_ref_nw4,poly_ref_nw5,
237##                poly_ref_nw6,poly_ref_nw7,poly_ref_nw8,poly_all],'poly_pic')
238
239interior_regions = [[poly_dampier,res_dampier], 
240                    [poly_karratha,res_karratha],[poly_karratha_town,res_karratha_town],
241                    [poly_delambre,res_delambre],[poly_coast,res_coast],
242                    [poly_facility,res_facility],
243                    #[poly_NWislands,res_NWislands],
244                    [poly_island0,res_island0],[poly_island1,res_island0],
245                    [poly_island2,res_island0],[poly_island3,res_island0],
246                    [poly_island4,res_island0],[poly_island5,res_island0],
247                    [poly_island6,res_island0],[poly_island7,res_island0],
248                    [poly_island8,res_island0],[poly_ref_nw4,res_ref_nw4],
249                    [poly_ref_nw5,res_ref_nw5],[poly_ref_nw6,res_ref_nw6],
250                    [poly_ref_nw7,res_ref_nw7],[poly_ref_nw8,res_ref_nw8]]
251
252interior_regions_test = [[poly_dampier,res_dampier], 
253                    [poly_karratha,res_karratha],[poly_karratha_town,res_karratha_town]]
254
255trigs_min = number_mesh_triangles(interior_regions_test, poly_all, res_poly_all)
256
257print 'min number triangles', trigs_min
258
259###################################################################
260# Clipping regions for export to asc and regions for clipping data
261###################################################################
262
263poly_bathy = read_polygon(polygons_dir+'mainland_only.csv')
264
265# exporting asc grid - Dampier
266e_min_area = 474000
267e_max_area = 480000
268n_min_area = 7719000
269n_max_area = 7725000
270
271# exporting asc grid - Karratha
272##e_min_area =
273##e_max_area =
274##n_min_area =
275##n_max_area =
276
277"""
278# used in the CIPMA 2006 scenario
279# CIPMA point of interest
280cipma_latitude = -20.588456
281cipma_longitude = 116.771527
282
283k0 = [cipma_latitude-0.02, cipma_longitude-0.02]
284k1 = [cipma_latitude-0.02, cipma_longitude+0.02]
285k2 = [cipma_latitude+0.02, cipma_longitude+0.02]
286k3 = [cipma_latitude+0.02, cipma_longitude-0.02]
287
288cipma_polygon, zone = convert_from_latlon_to_utm([k0, k1, k2, k3])
289assert zone == refzone
290
291poly_facility = read_polygon(polygons_dir+'facility.csv')
292poly_pipeline = read_polygon(polygons_dir+'pipeline2.csv')
293poly_interior = read_polygon(polygons_dir+'interior.csv')
294poly_coast = read_polygon(polygons_dir+'coast_final.csv')
295clip_poly_e = read_polygon(polygons_dir+'gap_e.csv')
296clip_poly_nw = read_polygon(polygons_dir+'gap_nw.csv')
297clip_poly_mid_w = read_polygon(polygons_dir+'gap_mid_w.cvs')
298clip_poly_mid_e = read_polygon(polygons_dir+'gap_mid_e.cvs')
299
300# exporting asc grid
301e_min_area = 474000
302e_max_area = 480000
303n_min_area = 7719000
304n_max_area = 7725000
305
306# used in the original 2005 scenario
307#Interior regions
308karratha_south = degminsec2decimal_degrees(-20,44,0)
309karratha_north = degminsec2decimal_degrees(-20,42,0)
310karratha_west = degminsec2decimal_degrees(116,48,0)
311karratha_east = degminsec2decimal_degrees(116,53,30)
312
313k0 = [karratha_south, karratha_west]
314k1 = [karratha_south, karratha_east]
315k2 = [karratha_north, karratha_east]
316k3 = [karratha_north, karratha_west]   
317
318karratha_polygon, zone = convert_from_latlon_to_utm([k0, k1, k2, k3])
319assert zone == refzone
320
321#Interior regions
322dampier_south = degminsec2decimal_degrees(-20,40,0)
323dampier_north = degminsec2decimal_degrees(-20,38,10)
324dampier_west = degminsec2decimal_degrees(116,43,0)
325dampier_east = degminsec2decimal_degrees(116,45,0)
326
327d0 = [dampier_south, dampier_west]
328d1 = [dampier_south, dampier_east]
329d2 = [dampier_north, dampier_east]
330d3 = [dampier_north, dampier_west]   
331
332dampier_polygon, zone = convert_from_latlon_to_utm([d0, d1, d2, d3])
333assert zone == refzone
334
335#Interior regions
336refinery_south = degminsec2decimal_degrees(-20,37,50)
337refinery_north = degminsec2decimal_degrees(-20,36,0)
338refinery_west = degminsec2decimal_degrees(116,44,0)
339refinery_east = degminsec2decimal_degrees(116,46,10)
340
341d0 = [refinery_south, refinery_west]
342d1 = [refinery_south, refinery_east]
343d2 = [refinery_north, refinery_east]
344d3 = [refinery_north, refinery_west]   
345
346refinery_polygon, zone = convert_from_latlon_to_utm([d0, d1, d2, d3])
347assert zone == refzone
348
349#Interior region around 468899, 7715177:
350#lat (-20, 39, 44.93753), lon (116, 42, 5.09106)
351
352point_south = degminsec2decimal_degrees(-20,39,46)
353point_north = degminsec2decimal_degrees(-20,39,42)
354point_west = degminsec2decimal_degrees(116,42,0)
355point_east = degminsec2decimal_degrees(116,42,10)
356
357d0 = [point_south, point_west]
358d1 = [point_south, point_east]
359d2 = [point_north, point_east]
360d3 = [point_north, point_west]   
361
362point_polygon, zone = convert_from_latlon_to_utm([d0, d1, d2, d3])
363assert zone == refzone
364
365#Neils areas around interesting points
366neil1_point1 = [degminsec2decimal_degrees(-20,35,34),
367                degminsec2decimal_degrees(116,45,18)]
368neil1_point2 = [degminsec2decimal_degrees(-20,36,15),
369                degminsec2decimal_degrees(116,46,18)]
370neil1_point3 = [degminsec2decimal_degrees(-20,35,9),
371                degminsec2decimal_degrees(116,47,17)]
372neil1_point4 = [degminsec2decimal_degrees(-20,34,26),
373                degminsec2decimal_degrees(116,46,17)]
374
375neil1_polygon, zone = convert_from_latlon_to_utm([neil1_point1,
376                                                  neil1_point2,
377                                                  neil1_point3,
378                                                  neil1_point4])
379assert zone == refzone
380
381neil2_point1 = [degminsec2decimal_degrees(-20,39,36),
382                degminsec2decimal_degrees(116,41,33)]
383neil2_point2 = [degminsec2decimal_degrees(-20,40,10),
384                degminsec2decimal_degrees(116,42,13)]
385neil2_point3 = [degminsec2decimal_degrees(-20,38,39),
386                degminsec2decimal_degrees(116,43,49)]
387neil2_point4 = [degminsec2decimal_degrees(-20,38,5),
388                degminsec2decimal_degrees(116,43,9)]
389
390neil2_polygon, zone = convert_from_latlon_to_utm([neil2_point1,
391                                                  neil2_point2,
392                                                  neil2_point3,
393                                                  neil2_point4])
394assert zone == refzone
395
396#Withnell bay
397wb_point1 = [degminsec2decimal_degrees(-20,35,34),
398                degminsec2decimal_degrees(116,45,18)]
399wb_point2 = [degminsec2decimal_degrees(-20,36,15),
400                degminsec2decimal_degrees(116,46,18)]
401wb_point3 = [degminsec2decimal_degrees(-20,35,9),
402                degminsec2decimal_degrees(116,47,17)]
403wb_point4 = [degminsec2decimal_degrees(-20,34,26),
404                degminsec2decimal_degrees(116,46,17)]
405
406wb_polygon, zone = convert_from_latlon_to_utm([wb_point1, wb_point2,
407                                               wb_point3, wb_point4])
408assert zone == refzone
409
410#Larger Withnell bay
411lwb_point1 = [degminsec2decimal_degrees(-20,35,59),
412                degminsec2decimal_degrees(116,42,00)]
413lwb_point2 = [degminsec2decimal_degrees(-20,36,50),
414                degminsec2decimal_degrees(116,46,50)]
415lwb_point3 = [degminsec2decimal_degrees(-20,34,00),
416                degminsec2decimal_degrees(116,47,39)]
417lwb_point4 = [degminsec2decimal_degrees(-20,33,00),
418                degminsec2decimal_degrees(116,42,50)]
419
420lwb_polygon, zone = convert_from_latlon_to_utm([lwb_point1, lwb_point2,
421                                                lwb_point3, lwb_point4])
422                                                     
423assert zone == refzone
424"""
Note: See TracBrowser for help on using the repository browser.