source: production/sydney_2006/project.py @ 3535

Last change on this file since 3535 was 3535, checked in by duncan, 18 years ago

change imports so reflect the new structure

File size: 16.3 KB
Line 
1"""Common filenames and locations for topographic data, meshes and outputs.
2Also includes origin for slump scenario.
3"""
4
5from os import sep, environ
6from os.path import expanduser
7from anuga.utilities.polygon import read_polygon
8import sys
9
10from anuga.pmesh.create_mesh import convert_points_from_latlon_to_utm
11               
12
13#Making assumptions about the location of scenario data
14state = 'new_south_wales'
15scenario_dir_name = 'sydney_tsunami_scenario_2006'
16# revised 100m data
17coarsename = 'bathyland100' # get from Neil/Ingo (DEM or topo data)
18# revised 25m data
19finename = 'bathy_dem25' # get from Neil/Ingo (DEM or topo data) Wed 25 Jan
20
21# creating easting and northing max and min for fine data - region of interest
22eastingmin = 332090
23eastingmax = 347500
24northingmin = 6246250
25northingmax = 6264100
26
27# creating easting and northing max and min for export viz purposes
28eminviz = 318000
29emaxviz = 351000
30nminviz = 6231000
31nmaxviz = 6283000
32
33basename = 'slump_duncan'
34basename2 = 'slump_1000res_1min'
35basename4 = 'slump_poly_ingo_test'
36
37if sys.platform == 'win32':
38    home = environ['INUNDATIONHOME']     #Sandpit's parent dir
39else:   
40    home = expanduser('~')
41   
42#Derive subdirectories and filenames
43meshdir = home+sep+state+scenario_dir_name+sep+'anuga'+sep+'meshes'+sep
44datadir = home+sep+state+scenario_dir_name+sep+'anuga'+sep+'topographies'+sep
45outputdir = home+sep+state+scenario_dir_name+sep+'anuga'+sep+'output'+sep
46polygondir = home+sep+state+scenario_dir_name+sep+'anuga'+sep+'polygons'+sep
47gaugedir = home+sep+state+scenario_dir_name+sep+'anuga'+sep+'gauges'+sep
48#reportdir = home+sep+state+scenario_dir_name+sep+'reports'+sep
49
50meshname = meshdir + basename
51meshelevname = meshdir + 'test_elev.tsh'
52meshname4 = meshdir + basename4
53coarsedemname = datadir + coarsename
54finedemname = datadir + finename
55combineddemname = datadir + 'sydneytopo'
56outputname = outputdir + basename  #Used by post processing
57outputname2 = outputdir+sep+'Coast_Polygons'+sep+'res1000'+sep+basename2  #Used by post processing
58
59#csv file of coastline 50m epsilon belt
60#manly_polygonname = polygondir + 'manly_polygon_UTM56_coarse'
61#manly_polygon = read_polygon(manly_polygonname + '.csv')
62#print manly_polygon
63
64gauge_filename = gaugedir + 'sydney_gauges.xya'
65buildings_filename = gaugedir + 'all_bld_ind.csv'
66#gauge_filename = gaugedir + 'sydney_gauges_test.xya'
67gauge_outname = gaugedir + 'gauges_max_output.xya'
68#gauge_filename = gaugedir + 'nest_gauges_Manly.xya'
69#gauge_filename = gaugedir + 'west_of_quay_yprofile.xya'
70#gauge_filename = gaugedir + 'GA_gauge.csv' # from Benfield
71#gauge_outname = gaugedir + 'gauges_max_output_next.xya'
72gaugetimeseries = gaugedir + 'gauges_time_series'
73polygonptsfile = polygondir + 'poly'
74integraltimeseries = outputdir + 'integral_time_series'
75integraltimeseries2 = outputdir + 'integral_time_series_move_origin'
76
77
78#Georeferencing
79from anuga.coordinate_transforms.redfearn import degminsec2decimal_degrees
80
81# define clipping polygon
82south = degminsec2decimal_degrees(-34,05,0)
83north = degminsec2decimal_degrees(-33,33,0)
84west = degminsec2decimal_degrees(151,1,0)
85east = degminsec2decimal_degrees(151,30,0)
86p0 = [south, west]
87p1 = [south, east]
88p2 = [north, east]
89p3 = [north, west]
90   
91polygonall, zone = convert_points_from_latlon_to_utm([p0, p1, p2, p3])
92refzone = zone
93
94print 'Got refzone', refzone
95
96# original clipping
97dsouth = degminsec2decimal_degrees(-34,05,0)
98dsouth2 = degminsec2decimal_degrees(-34,01,0)
99dnorth = degminsec2decimal_degrees(-33,33,0)
100dnorth1 = degminsec2decimal_degrees(-33,40,0)
101dnorth2 = degminsec2decimal_degrees(-33,58,30)
102dnorth3 = degminsec2decimal_degrees(-33,46,0)
103dwest = degminsec2decimal_degrees(151,1,0)
104deast1 = degminsec2decimal_degrees(151,20,0)
105deast2 = degminsec2decimal_degrees(151,48,0)
106deast3 = degminsec2decimal_degrees(151,10,0)
107deast4 = degminsec2decimal_degrees(151,9,0)
108
109dp0 = [dsouth, dwest]
110dp1 = [dsouth, deast1]
111dp2 = [dnorth2, deast2]
112dp3 = [dnorth1, deast2]
113dp4 = [dnorth, deast1]
114dp5 = [dnorth, deast4]
115dp6 = [dnorth3, deast3]
116dp7 = [dnorth3, dwest]
117dp8 = [dnorth, dwest]
118
119dp12 = [dsouth2, deast2]
120dp34 = [dnorth, deast2]
121diffpolygonall, zone = convert_points_from_latlon_to_utm([dp0, dp1, dp2, dp3, dp4, dp5, dp6, dp7])
122# used for new tests 4 April 2006 (ensure slump contained in domain)
123diffpolygonall2, zone = convert_points_from_latlon_to_utm([dp0, dp1, dp12, dp34, dp4, dp5, dp6, dp7])
124
125# clipping used for look at increasingly finer resolution
126j0 = [338000, 6243000]
127j1 = [365000, 6243000]
128j2 = [365000, 6273000]
129j3 = [338000, 6273000]
130   
131diffpolygonall_test = [j0, j1, j2, j3]
132# clipping used for further investigation of black screen of death
133j0 = [328000, 6255000]
134j1 = [355000, 6255000]
135j2 = [355000, 6270000]
136j3 = [328000, 6270000]
137
138diffpolygonall_test2 = [j0, j1, j2, j3]
139
140# clipping used for demo purposes to fit around poly3 and poly4 (Ingo's files)
141#j0 = [318000, 6249000]
142#j1 = [387000, 6249000]
143#j2 = [387000, 6270000]
144#j3 = [318000, 6270000]
145
146#demopoly = [j0, j1, j2, j3]
147
148# second demo poly which isn't rectangular
149j0 = [385000, 6280000]
150j1 = [360000, 6272500]
151j2 = [335000, 6272500]
152j3 = [330000, 6265000]
153j31 = [325000, 6260000]
154j4 = [316000, 6260000]
155j5 = [316000, 6247000]
156j6 = [350000, 6247000]
157j7 = [385000, 6238000]
158
159demopoly = [j0, j1, j2, j3, j31, j4, j5, j6, j7]
160
161# to put chunk back in
162#diffpolygonall = [dp0, dp1, dp2, dp3, dp4, dp8]
163
164# testing new interior regions 15 Feb 06
165# these worked OK it seemed, no warnings and resulting mesh looked fine.
166pp0 = [343965, 6273229]
167pp1 = [342984, 6270664]
168pp2 = [343950, 6270005]
169pp3 = [343853, 6270399]
170pp4 = [343383, 6270925]
171pp5 = [343756, 6271926]
172pp6 = [344037, 6272401]
173pp7 = [344405, 6272503]
174pp8 = [344226, 6273244]
175
176poly1 = [pp0, pp1, pp2, pp3, pp4, pp5, pp6, pp7, pp8]
177
178qp0 = [343494.0, 6270650.0]
179qp1 = [343337.0, 6270303.0]
180qp2 = [343466.0, 6270228.0]
181qp3 = [343139.0, 6269901.0]
182qp4 = [342472.0, 6268573.0]
183qp5 = [342111.0, 6267115.0]
184qp6 = [342479.0, 6266121.0]
185qp7 = [342860.0, 6266386.0]
186qp8 = [342635.0, 6267438.0]
187qp9 = [343105.0, 6269070.0]
188qp10 = [343548.0, 6269567.0]
189qp11 = [343487.0, 6269928.0]
190qp12 = [343991.0, 6270269.0]
191
192poly2 = [qp0, qp1, qp2, qp3, qp4, qp5, qp6, qp7, qp8, qp9, qp10, qp11, qp12]
193#poly2 = [qp5, qp6, qp7, qp8, qp9, qp10]
194
195# test from Jane drawing
196# north
197np1 = [318200, 6253000]
198np2 = [327000, 6250500]
199np2east = [327200, 6250500]
200np2west = [326800, 6250500]
201np3 = [333000, 6249000]
202np4 = [342000, 6249000]
203np5 = [342000, 6255000]
204np6 = [343000, 6256000]
205np7 = [344000, 6258000]
206np8 = [343000, 6260300]
207#np8 = [343000, 6260000]
208np9 = [343000, 6264000]
209np10 = [345000, 6265000]
210np10_2 = [344500, 6265000]
211np10_3 = [344200, 6265000]
212np11 = [343000, 6266000]
213np12 = [344000, 6269000]
214np13 = [344000, 6272000]
215np14 = [342000, 6272000]
216np15 = [339000, 6269000]
217np16 = [339000, 6264000]
218#np17 = [332000, 6264000]
219np17 = [332500, 6262000]
220np18 = [334000, 6254000]
221np19 = [329000, 6257000]
222np20 = [330000, 6259000]
223#np21 = [327000, 6262000]
224np21 = [327000, 6259000]
225np22 = [327000, 6257000]
226np23 = [318200, 6257000]
227np24 = [335000, 6261000]
228np25 = [336000, 6262000]
229np26 = [338000, 6262000]
230np27 = [340000, 6264000]
231np28 = [334000, 6250000]
232np29 = [336000, 6250000]
233np30 = [341500, 6250000]
234np31 = [327000, 6254000]
235np32 = [323000, 6257000]
236np33 = [335000, 6255000]
237np34 = [337000, 6256000]
238np35 = [338000, 6255000]
239np36 = [320000, 6254000]
240np37 = [322000, 6252000]
241np38 = [324000, 6253000]
242np39 = [325000, 6251000]
243nptest = [339000, 6266000]
244nptest2 = [344000, 6265000]
245
246#testpoly1 = [np1, \
247#             np36, np37, np38, np39, \
248#             np2, np3, \
249#             np28, np29, np30, \
250#             np5, np6, np7, np8, np9, np10, \
251#             np11, np12, np13, np14, np15, \
252#             np27, np26, np25, np24, \
253#             np17, \
254#             np33, np34, np35, \
255#             np18, np19, np20, \
256#             np21, np22, \
257#             np31, np32, \
258#             np23]
259
260# reduce this polygon for black screen of death testing 24/04/06
261#newpoly1_refine = [nptest, np27, np16] # worked
262#newpoly1_refine = [np9, np10, np11] # worked
263#newpoly1_refine = [np9, np10, np11, nptest2] # worked
264#newpoly1_refine = [np9, np10_2, np11, nptest2] # worked
265newpoly1_refine = [np9, np10_3, np11, nptest2] # worked
266
267# used for refined interior regions and fine mesh
268newpoly1 = [np28, np29, np30, \
269            np5, np6, np7, np8, np9, np10, \
270            np11, np12, np13, np14, np15, \
271            np27, np26, np25, np24, np17, \
272            np33, \
273            np19, np20, np21, np22, np31, np32, \
274            np36, np37, np38, np39, np2, np3]
275
276# last two lines for second run
277
278north1 = [np2, np3, np30, np5, np24, np17, np33, np19, np20, np21, \
279          np31, np32, np37, np38, np39]
280north2 = [np8, np9, np19, np11, np12, np13, np14, np15, np27, np27, np25]
281npinsert = [332000, 6255000]
282parrariver = [np3, np18, npinsert, np19, np32, np36, np38, np39, np2] #first run
283#parrariver = [np32, np36, np37, np38, np39, np2west, np31] #second run
284
285# south
286sp1 = [ 328000, 6231000]
287sp2 = [335000, 6231000]
288sp3 = [338000, 6235000]
289sp4 = [340000, 6240000]
290sp5 = [340000, 6244000]
291sp6 = [342000, 6248800]
292sp7 = [340000, 6248800]
293sp8 = [338000, 6244000]
294sp9 = [338000, 6237000]
295sp10 = [334000, 6242000]
296sp11 = [331000, 6245000]
297sp12 = [326000, 6247000]
298sp13 = [325000, 6246000]
299sp14 = [329000, 6243000]
300sp15 = [328000, 6237000]
301sp16 = [337000, 6236000]
302sp17 = [330000, 6237000]
303sp18 = [330000, 6239000]
304sp19 = [332000, 6239000]
305sp20 = [334000, 6240000]
306sp21 = [334000, 6238000]
307sp22 = [337000, 6236500]
308sp23 = [339000, 6236000]
309
310#original
311south1 = [sp1, sp2, sp3, \
312             sp16, sp17, sp18, sp19, sp20, sp21, sp22, sp23, \
313             sp4, sp5, sp6, sp7, sp8, sp9, sp10, \
314             sp11, sp12, sp13, sp14, sp15] 
315
316m1 = [340000, 6256000]
317m2 = [342800, 6256000]
318m3 = [342800, 6261000]
319m4 = [340000, 6260000]
320   
321finepolymanly = [m1, m2, m3, m4]
322
323q1 = [333000, 6250000]
324q2 = [340000, 6250000]
325q3 = [340000, 6254000]
326q4 = [333000, 6255000]
327   
328finepolyquay = [q1, q2, q3, q4]
329
330# refined 13/03/06 to look at effect of finer resolution on no IC scenario
331q2south = [340000, 6249000]
332newpoly1 = [np28, np29, q2south, sp7, sp8, sp5, sp6, np30, \
333            np5, np6, np7, np8, np9, np10, \
334            np11, np12, np13, np14, np15, \
335            np27, np26, np25, np24, np17, \
336            np33, \
337            np19, np20, np21, np22, np31, np32, \
338            np36, np37, np38, np39, np2, np3]
339
340#Interior regions - the Harbour
341harbour_1x = degminsec2decimal_degrees(-33,51,0)
342harbour_1y = degminsec2decimal_degrees(151,2,30)
343harbour_12x = degminsec2decimal_degrees(-33,51,0)
344harbour_12y = degminsec2decimal_degrees(151,5,0)
345harbour_13x = degminsec2decimal_degrees(-33,52,15)
346harbour_13y = degminsec2decimal_degrees(151,5,0)
347harbour_2x = degminsec2decimal_degrees(-33,53,0)
348harbour_2y = degminsec2decimal_degrees(151,17,20)
349harbour_3x = degminsec2decimal_degrees(-33,47,0)
350harbour_3y = degminsec2decimal_degrees(151,20,30)
351harbour_4x = degminsec2decimal_degrees(-33,47,50)
352harbour_4y = degminsec2decimal_degrees(151,8,10)
353harbour_5x = degminsec2decimal_degrees(-33,48,10)
354harbour_5y = degminsec2decimal_degrees(151,8,0)
355harbour_6x = degminsec2decimal_degrees(-33,49,0)
356harbour_6y = degminsec2decimal_degrees(151,2,30)
357harbour_7x = degminsec2decimal_degrees(-33,34,30)
358harbour_7y = degminsec2decimal_degrees(151,20,20)
359harbour_8x = degminsec2decimal_degrees(-33,33,30)
360harbour_8y = degminsec2decimal_degrees(151,17,0)
361harbour_9x = degminsec2decimal_degrees(-33,45,30)
362harbour_9y = degminsec2decimal_degrees(151,17,0)
363harbour_10x = degminsec2decimal_degrees(-33,45,10)
364harbour_10y = degminsec2decimal_degrees(151,11,40)
365harbour_11x = degminsec2decimal_degrees(-33,45,10)
366harbour_11y = degminsec2decimal_degrees(151,11,40)
367harbour_14x = degminsec2decimal_degrees(-33,49,10)
368harbour_14y = degminsec2decimal_degrees(151,11,40)
369harbour_15x = degminsec2decimal_degrees(-33,48,55)
370harbour_15y = degminsec2decimal_degrees(151,2,30)
371
372k02 = [harbour_1x, harbour_1y]
373k12 = [harbour_2x, harbour_2y]
374k22 = [harbour_3x, harbour_3y]
375k32 = [harbour_4x, harbour_4y]
376k42 = [harbour_5x, harbour_5y]
377k52 = [harbour_6x, harbour_6y]
378k62 = [harbour_7x, harbour_7y]
379k72 = [harbour_8x, harbour_8y]
380k82 = [harbour_9x, harbour_9y]
381k92 = [harbour_10x, harbour_10y]
382k102 = [harbour_11x, harbour_11y]
383k112 = [harbour_12x, harbour_12y]
384k122 = [harbour_13x, harbour_13y]
385k132 = [harbour_14x, harbour_14y]
386k142 = [harbour_15x, harbour_15y]
387
388harbour_polygon_2, zone = convert_points_from_latlon_to_utm([k02, k112, k122, k12, k22, k62, k72, k82, k102, k42, k52]) #worked
389assert zone == refzone
390
391
392#Interior region - Botany Bay
393bb_1x = degminsec2decimal_degrees(-34,3,0)
394bb_1y = degminsec2decimal_degrees(151,2,30)
395bb_10x = degminsec2decimal_degrees(-34,3,0)
396bb_10y = degminsec2decimal_degrees(151,8,0)
397bb_2x = degminsec2decimal_degrees(-34,3,0)
398bb_2y = degminsec2decimal_degrees(151,14,0)
399bb_3x = degminsec2decimal_degrees(-33,53,30)
400bb_3y = degminsec2decimal_degrees(151,17,20)
401bb_4x = degminsec2decimal_degrees(-33,53,0)
402bb_4y = degminsec2decimal_degrees(151,8,0)
403bb_5x = degminsec2decimal_degrees(-33,57,30)
404bb_5y = degminsec2decimal_degrees(151,8,0)
405bb_6x = degminsec2decimal_degrees(-33,57,30)
406bb_6y = degminsec2decimal_degrees(151,2,30) 
407bb_7x = degminsec2decimal_degrees(-33,53,30)
408bb_7y = degminsec2decimal_degrees(151,12,30) 
409bb_8x = degminsec2decimal_degrees(-33,55,20)
410bb_8y = degminsec2decimal_degrees(151,8,0) 
411bb_9x = degminsec2decimal_degrees(-33,55,20)
412bb_9y = degminsec2decimal_degrees(151,12,30) 
413
414j02 = [bb_1x, bb_1y]
415j12 = [bb_2x, bb_2y]
416j22 = [bb_3x, bb_3y]
417j32 = [bb_4x, bb_4y]
418j42 = [bb_5x, bb_5y]
419j52 = [bb_6x, bb_6y]
420j62 = [bb_7x, bb_7y]
421j72 = [bb_8x, bb_8y]
422j82 = [bb_9x, bb_9y]
423j92 = [bb_10x, bb_10y]
424
425botanybay_polygon_2, zone = convert_points_from_latlon_to_utm([j92, j12, j22, j62, j82, j72, j42]) # worked
426
427
428# close to botany bay opening (340000,6236000)
429# x0 = 25964
430# y0 = 11049
431# around 10km from botany bay opening (350000,6236000)
432# x0 = 35964
433# y0 = 11049
434# around 21km from botany bay opening (361000,6236000)
435#x0 = 46964
436#y0 = 11049
437
438# not used for sydney scenario, original interior regions listed though
439# setting up problem area for doing just around the harbour
440hsouth = degminsec2decimal_degrees(-33,54,0)
441hnorth = degminsec2decimal_degrees(-33,48,0)
442hwest = degminsec2decimal_degrees(151,0,0)
443heast = degminsec2decimal_degrees(151,30,0)
444
445hp0 = [hsouth, hwest]
446hp1 = [hsouth, heast]
447hp2 = [hnorth, heast]
448hp3 = [hnorth, hwest]
449polygon_h, zone = convert_points_from_latlon_to_utm([hp0, hp1, hp2, hp3])
450
451#Interior regions - the Harbour - take 1
452harbour_south = degminsec2decimal_degrees(-33,53,0)
453harbour_north = degminsec2decimal_degrees(-33,47,0)
454harbour_west = degminsec2decimal_degrees(151,5,0)
455harbour_east = degminsec2decimal_degrees(151,19,0)
456
457#harbour_south1 = degminsec2decimal_degrees(-33,53,0)
458#harbour_south2 = degminsec2decimal_degrees(-33,52,0)
459#harbour_north1 = degminsec2decimal_degrees(-33,45,0)
460#harbour_north2 = degminsec2decimal_degrees(-33,48,0)
461#harbour_west = degminsec2decimal_degrees(151,5,0)
462#harbour_east = degminsec2decimal_degrees(151,19,0)
463
464k0 = [harbour_south, harbour_west]
465k1 = [harbour_south, harbour_east]
466k2 = [harbour_north, harbour_east]
467k3 = [harbour_north, harbour_west]   
468
469harbour_polygon, zone = convert_points_from_latlon_to_utm([k0, k1, k2, k3])
470
471# setting up problem area for doing just around Botany Bay
472bsouth = degminsec2decimal_degrees(-33,56,0)
473bnorth = degminsec2decimal_degrees(-34,3,0)
474bwest = degminsec2decimal_degrees(151,0,0)
475beast = degminsec2decimal_degrees(151,30,0)
476
477bp0 = [bsouth, bwest]
478bp1 = [bsouth, beast]
479bp2 = [bnorth, beast]
480bp3 = [bnorth, bwest]
481polygon_bb, zone = convert_points_from_latlon_to_utm([bp0, bp1, bp2, bp3])
482
483#Interior region - Botany Bay - take 1
484botanybay_south = degminsec2decimal_degrees(-33,58,0)
485botanybay_north = degminsec2decimal_degrees(-34,1,0)
486botanybay_west = degminsec2decimal_degrees(151,5,0)
487botanybay_east = degminsec2decimal_degrees(151,18,0)
488
489j0 = [botanybay_south, botanybay_west]
490j1 = [botanybay_south, botanybay_east]
491j2 = [botanybay_north, botanybay_east]
492j3 = [botanybay_north, botanybay_west]   
493
494botanybay_polygon, zone = convert_points_from_latlon_to_utm([j0, j1, j2, j3])
495assert zone == refzone
496
497#x0 = 28964 + 42000
498#y0 = 30049
499#slump_origin = [x0+314036.58727982, y0+6224951.2960092] #Absolute UTM
500slump_origin = [385000.0, 6255000.0] #Absolute UTM
501# move 10km west
502slump_origin2 = [375000.0, 6255000.0] #Absolute UTM
503
504a = [340000, 6255000]
505b = [340000, 6270000]
506c = [318000, 6255000]
507d = [318000, 6270000]
508e = [316000, 6255000]
509f = [395000, 6255000]
510g = [355000, 6280000]
511h = [355000, 6224000]
512test_pts = [a, b, c, d, e, f, g, h]
513test_elev = [1.0, 4.0, 3.0, 0.1, 5, -10.0, -20, -15]
Note: See TracBrowser for help on using the repository browser.