source: production/sydney_2006/project.py @ 3514

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

Hi all,
I'm doing a change in the anuga structure, moving the code to

\anuga_core\source\anuga

After you have done an svn update, the PYTHONPATH has to be changed to;
PYTHONPATH = anuga_core/source/

This is part of changes required to make installation of anuga quicker and reducing the size of our sandpits.

If any imports are broken, try fixing them. With adding anuga. to them for example. If this seems to have really broken things, email/phone me.

Cheers
Duncan

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 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 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.