"""Common filenames and locations for topographic data, meshes and outputs. Also includes origin for slump scenario. """ from os import sep, environ from os.path import expanduser from anuga.utilities.polygon import read_polygon import sys from anuga.coordinate_transforms.redfearn import degminsec2decimal_degrees, convert_from_latlon_to_utm #Making assumptions about the location of scenario data state = 'new_south_wales' scenario_dir_name = 'sydney_tsunami_scenario_2006' # revised 100m data coarsename = 'bathyland100' # get from Neil/Ingo (DEM or topo data) # revised 25m data finename = 'bathy_dem25' # get from Neil/Ingo (DEM or topo data) Wed 25 Jan # creating easting and northing max and min for fine data - region of interest eastingmin = 332090 eastingmax = 347500 northingmin = 6246250 northingmax = 6264100 # creating easting and northing max and min for export viz purposes eminviz = 318000 emaxviz = 351000 nminviz = 6231000 nmaxviz = 6283000 basename = 'slump_duncan' basename2 = 'slump_1000res_1min' basename4 = 'slump_poly_ingo_test' if sys.platform == 'win32': home = environ['INUNDATIONHOME'] #Sandpit's parent dir else: home = expanduser('~') # INUNDATIONHOME is the inundation directory, not the data directory. home += sep +'data' #Derive subdirectories and filenames meshdir = home+sep+state+scenario_dir_name+sep+'anuga'+sep+'meshes'+sep datadir = home+sep+state+scenario_dir_name+sep+'anuga'+sep+'topographies'+sep outputdir = home+sep+state+scenario_dir_name+sep+'anuga'+sep+'output'+sep polygondir = home+sep+state+scenario_dir_name+sep+'anuga'+sep+'polygons'+sep gaugedir = home+sep+state+scenario_dir_name+sep+'anuga'+sep+'gauges'+sep #reportdir = home+sep+state+scenario_dir_name+sep+'reports'+sep meshname = meshdir + basename meshelevname = meshdir + 'test_elev.tsh' meshname4 = meshdir + basename4 coarsedemname = datadir + coarsename finedemname = datadir + finename combineddemname = datadir + 'sydneytopo' outputname = outputdir + basename #Used by post processing outputname2 = outputdir+sep+'Coast_Polygons'+sep+'res1000'+sep+basename2 #Used by post processing #csv file of coastline 50m epsilon belt #manly_polygonname = polygondir + 'manly_polygon_UTM56_coarse' #manly_polygon = read_polygon(manly_polygonname + '.csv') #print manly_polygon gauge_filename = gaugedir + 'sydney_gauges.xya' buildings_filename = gaugedir + 'all_bld_ind.csv' #gauge_filename = gaugedir + 'sydney_gauges_test.xya' gauge_outname = gaugedir + 'gauges_max_output.xya' #gauge_filename = gaugedir + 'nest_gauges_Manly.xya' #gauge_filename = gaugedir + 'west_of_quay_yprofile.xya' #gauge_filename = gaugedir + 'GA_gauge.csv' # from Benfield #gauge_outname = gaugedir + 'gauges_max_output_next.xya' gaugetimeseries = gaugedir + 'gauges_time_series' polygonptsfile = polygondir + 'poly' integraltimeseries = outputdir + 'integral_time_series' integraltimeseries2 = outputdir + 'integral_time_series_move_origin' #Georeferencing from anuga.coordinate_transforms.redfearn import degminsec2decimal_degrees # define clipping polygon south = degminsec2decimal_degrees(-34,05,0) north = degminsec2decimal_degrees(-33,33,0) west = degminsec2decimal_degrees(151,1,0) east = degminsec2decimal_degrees(151,30,0) p0 = [south, west] p1 = [south, east] p2 = [north, east] p3 = [north, west] polygonall, zone = convert_from_latlon_to_utm([p0, p1, p2, p3]) refzone = zone print 'Got refzone', refzone # original clipping dsouth = degminsec2decimal_degrees(-34,05,0) dsouth2 = degminsec2decimal_degrees(-34,01,0) dnorth = degminsec2decimal_degrees(-33,33,0) dnorth1 = degminsec2decimal_degrees(-33,40,0) dnorth2 = degminsec2decimal_degrees(-33,58,30) dnorth3 = degminsec2decimal_degrees(-33,46,0) dwest = degminsec2decimal_degrees(151,1,0) deast1 = degminsec2decimal_degrees(151,20,0) deast2 = degminsec2decimal_degrees(151,48,0) deast3 = degminsec2decimal_degrees(151,10,0) deast4 = degminsec2decimal_degrees(151,9,0) dp0 = [dsouth, dwest] dp1 = [dsouth, deast1] dp2 = [dnorth2, deast2] dp3 = [dnorth1, deast2] dp4 = [dnorth, deast1] dp5 = [dnorth, deast4] dp6 = [dnorth3, deast3] dp7 = [dnorth3, dwest] dp8 = [dnorth, dwest] dp12 = [dsouth2, deast2] dp34 = [dnorth, deast2] diffpolygonall, zone = convert_from_latlon_to_utm([dp0, dp1, dp2, dp3, dp4, dp5, dp6, dp7]) # used for new tests 4 April 2006 (ensure slump contained in domain) diffpolygonall2, zone = convert_from_latlon_to_utm([dp0, dp1, dp12, dp34, dp4, dp5, dp6, dp7]) # clipping used for look at increasingly finer resolution j0 = [338000, 6243000] j1 = [365000, 6243000] j2 = [365000, 6273000] j3 = [338000, 6273000] diffpolygonall_test = [j0, j1, j2, j3] # clipping used for further investigation of black screen of death j0 = [328000, 6255000] j1 = [355000, 6255000] j2 = [355000, 6270000] j3 = [328000, 6270000] diffpolygonall_test2 = [j0, j1, j2, j3] # clipping used for demo purposes to fit around poly3 and poly4 (Ingo's files) #j0 = [318000, 6249000] #j1 = [387000, 6249000] #j2 = [387000, 6270000] #j3 = [318000, 6270000] #demopoly = [j0, j1, j2, j3] # second demo poly which isn't rectangular j0 = [385000, 6280000] j1 = [360000, 6272500] j2 = [335000, 6272500] j3 = [330000, 6265000] j31 = [325000, 6260000] j4 = [316000, 6260000] j5 = [316000, 6247000] j6 = [350000, 6247000] j7 = [385000, 6238000] demopoly = [j0, j1, j2, j3, j31, j4, j5, j6, j7] # to put chunk back in #diffpolygonall = [dp0, dp1, dp2, dp3, dp4, dp8] # testing new interior regions 15 Feb 06 # these worked OK it seemed, no warnings and resulting mesh looked fine. pp0 = [343965, 6273229] pp1 = [342984, 6270664] pp2 = [343950, 6270005] pp3 = [343853, 6270399] pp4 = [343383, 6270925] pp5 = [343756, 6271926] pp6 = [344037, 6272401] pp7 = [344405, 6272503] pp8 = [344226, 6273244] poly1 = [pp0, pp1, pp2, pp3, pp4, pp5, pp6, pp7, pp8] qp0 = [343494.0, 6270650.0] qp1 = [343337.0, 6270303.0] qp2 = [343466.0, 6270228.0] qp3 = [343139.0, 6269901.0] qp4 = [342472.0, 6268573.0] qp5 = [342111.0, 6267115.0] qp6 = [342479.0, 6266121.0] qp7 = [342860.0, 6266386.0] qp8 = [342635.0, 6267438.0] qp9 = [343105.0, 6269070.0] qp10 = [343548.0, 6269567.0] qp11 = [343487.0, 6269928.0] qp12 = [343991.0, 6270269.0] poly2 = [qp0, qp1, qp2, qp3, qp4, qp5, qp6, qp7, qp8, qp9, qp10, qp11, qp12] #poly2 = [qp5, qp6, qp7, qp8, qp9, qp10] # test from Jane drawing # north np1 = [318200, 6253000] np2 = [327000, 6250500] np2east = [327200, 6250500] np2west = [326800, 6250500] np3 = [333000, 6249000] np4 = [342000, 6249000] np5 = [342000, 6255000] np6 = [343000, 6256000] np7 = [344000, 6258000] np8 = [343000, 6260300] #np8 = [343000, 6260000] np9 = [343000, 6264000] np10 = [345000, 6265000] np10_2 = [344500, 6265000] np10_3 = [344200, 6265000] np11 = [343000, 6266000] np12 = [344000, 6269000] np13 = [344000, 6272000] np14 = [342000, 6272000] np15 = [339000, 6269000] np16 = [339000, 6264000] #np17 = [332000, 6264000] np17 = [332500, 6262000] np18 = [334000, 6254000] np19 = [329000, 6257000] np20 = [330000, 6259000] #np21 = [327000, 6262000] np21 = [327000, 6259000] np22 = [327000, 6257000] np23 = [318200, 6257000] np24 = [335000, 6261000] np25 = [336000, 6262000] np26 = [338000, 6262000] np27 = [340000, 6264000] np28 = [334000, 6250000] np29 = [336000, 6250000] np30 = [341500, 6250000] np31 = [327000, 6254000] np32 = [323000, 6257000] np33 = [335000, 6255000] np34 = [337000, 6256000] np35 = [338000, 6255000] np36 = [320000, 6254000] np37 = [322000, 6252000] np38 = [324000, 6253000] np39 = [325000, 6251000] nptest = [339000, 6266000] nptest2 = [344000, 6265000] #testpoly1 = [np1, \ # np36, np37, np38, np39, \ # np2, np3, \ # np28, np29, np30, \ # np5, np6, np7, np8, np9, np10, \ # np11, np12, np13, np14, np15, \ # np27, np26, np25, np24, \ # np17, \ # np33, np34, np35, \ # np18, np19, np20, \ # np21, np22, \ # np31, np32, \ # np23] # reduce this polygon for black screen of death testing 24/04/06 #newpoly1_refine = [nptest, np27, np16] # worked #newpoly1_refine = [np9, np10, np11] # worked #newpoly1_refine = [np9, np10, np11, nptest2] # worked #newpoly1_refine = [np9, np10_2, np11, nptest2] # worked newpoly1_refine = [np9, np10_3, np11, nptest2] # worked # used for refined interior regions and fine mesh newpoly1 = [np28, np29, np30, \ np5, np6, np7, np8, np9, np10, \ np11, np12, np13, np14, np15, \ np27, np26, np25, np24, np17, \ np33, \ np19, np20, np21, np22, np31, np32, \ np36, np37, np38, np39, np2, np3] # last two lines for second run north1 = [np2, np3, np30, np5, np24, np17, np33, np19, np20, np21, \ np31, np32, np37, np38, np39] north2 = [np8, np9, np19, np11, np12, np13, np14, np15, np27, np27, np25] npinsert = [332000, 6255000] parrariver = [np3, np18, npinsert, np19, np32, np36, np38, np39, np2] #first run #parrariver = [np32, np36, np37, np38, np39, np2west, np31] #second run # south sp1 = [ 328000, 6231000] sp2 = [335000, 6231000] sp3 = [338000, 6235000] sp4 = [340000, 6240000] sp5 = [340000, 6244000] sp6 = [342000, 6248800] sp7 = [340000, 6248800] sp8 = [338000, 6244000] sp9 = [338000, 6237000] sp10 = [334000, 6242000] sp11 = [331000, 6245000] sp12 = [326000, 6247000] sp13 = [325000, 6246000] sp14 = [329000, 6243000] sp15 = [328000, 6237000] sp16 = [337000, 6236000] sp17 = [330000, 6237000] sp18 = [330000, 6239000] sp19 = [332000, 6239000] sp20 = [334000, 6240000] sp21 = [334000, 6238000] sp22 = [337000, 6236500] sp23 = [339000, 6236000] #original south1 = [sp1, sp2, sp3, \ sp16, sp17, sp18, sp19, sp20, sp21, sp22, sp23, \ sp4, sp5, sp6, sp7, sp8, sp9, sp10, \ sp11, sp12, sp13, sp14, sp15] m1 = [340000, 6256000] m2 = [342800, 6256000] m3 = [342800, 6261000] m4 = [340000, 6260000] finepolymanly = [m1, m2, m3, m4] q1 = [333000, 6250000] q2 = [340000, 6250000] q3 = [340000, 6254000] q4 = [333000, 6255000] finepolyquay = [q1, q2, q3, q4] # refined 13/03/06 to look at effect of finer resolution on no IC scenario q2south = [340000, 6249000] newpoly1 = [np28, np29, q2south, sp7, sp8, sp5, sp6, np30, \ np5, np6, np7, np8, np9, np10, \ np11, np12, np13, np14, np15, \ np27, np26, np25, np24, np17, \ np33, \ np19, np20, np21, np22, np31, np32, \ np36, np37, np38, np39, np2, np3] #Interior regions - the Harbour harbour_1x = degminsec2decimal_degrees(-33,51,0) harbour_1y = degminsec2decimal_degrees(151,2,30) harbour_12x = degminsec2decimal_degrees(-33,51,0) harbour_12y = degminsec2decimal_degrees(151,5,0) harbour_13x = degminsec2decimal_degrees(-33,52,15) harbour_13y = degminsec2decimal_degrees(151,5,0) harbour_2x = degminsec2decimal_degrees(-33,53,0) harbour_2y = degminsec2decimal_degrees(151,17,20) harbour_3x = degminsec2decimal_degrees(-33,47,0) harbour_3y = degminsec2decimal_degrees(151,20,30) harbour_4x = degminsec2decimal_degrees(-33,47,50) harbour_4y = degminsec2decimal_degrees(151,8,10) harbour_5x = degminsec2decimal_degrees(-33,48,10) harbour_5y = degminsec2decimal_degrees(151,8,0) harbour_6x = degminsec2decimal_degrees(-33,49,0) harbour_6y = degminsec2decimal_degrees(151,2,30) harbour_7x = degminsec2decimal_degrees(-33,34,30) harbour_7y = degminsec2decimal_degrees(151,20,20) harbour_8x = degminsec2decimal_degrees(-33,33,30) harbour_8y = degminsec2decimal_degrees(151,17,0) harbour_9x = degminsec2decimal_degrees(-33,45,30) harbour_9y = degminsec2decimal_degrees(151,17,0) harbour_10x = degminsec2decimal_degrees(-33,45,10) harbour_10y = degminsec2decimal_degrees(151,11,40) harbour_11x = degminsec2decimal_degrees(-33,45,10) harbour_11y = degminsec2decimal_degrees(151,11,40) harbour_14x = degminsec2decimal_degrees(-33,49,10) harbour_14y = degminsec2decimal_degrees(151,11,40) harbour_15x = degminsec2decimal_degrees(-33,48,55) harbour_15y = degminsec2decimal_degrees(151,2,30) k02 = [harbour_1x, harbour_1y] k12 = [harbour_2x, harbour_2y] k22 = [harbour_3x, harbour_3y] k32 = [harbour_4x, harbour_4y] k42 = [harbour_5x, harbour_5y] k52 = [harbour_6x, harbour_6y] k62 = [harbour_7x, harbour_7y] k72 = [harbour_8x, harbour_8y] k82 = [harbour_9x, harbour_9y] k92 = [harbour_10x, harbour_10y] k102 = [harbour_11x, harbour_11y] k112 = [harbour_12x, harbour_12y] k122 = [harbour_13x, harbour_13y] k132 = [harbour_14x, harbour_14y] k142 = [harbour_15x, harbour_15y] harbour_polygon_2, zone = convert_from_latlon_to_utm([k02, k112, k122, k12, k22, k62, k72, k82, k102, k42, k52]) #worked assert zone == refzone #Interior region - Botany Bay bb_1x = degminsec2decimal_degrees(-34,3,0) bb_1y = degminsec2decimal_degrees(151,2,30) bb_10x = degminsec2decimal_degrees(-34,3,0) bb_10y = degminsec2decimal_degrees(151,8,0) bb_2x = degminsec2decimal_degrees(-34,3,0) bb_2y = degminsec2decimal_degrees(151,14,0) bb_3x = degminsec2decimal_degrees(-33,53,30) bb_3y = degminsec2decimal_degrees(151,17,20) bb_4x = degminsec2decimal_degrees(-33,53,0) bb_4y = degminsec2decimal_degrees(151,8,0) bb_5x = degminsec2decimal_degrees(-33,57,30) bb_5y = degminsec2decimal_degrees(151,8,0) bb_6x = degminsec2decimal_degrees(-33,57,30) bb_6y = degminsec2decimal_degrees(151,2,30) bb_7x = degminsec2decimal_degrees(-33,53,30) bb_7y = degminsec2decimal_degrees(151,12,30) bb_8x = degminsec2decimal_degrees(-33,55,20) bb_8y = degminsec2decimal_degrees(151,8,0) bb_9x = degminsec2decimal_degrees(-33,55,20) bb_9y = degminsec2decimal_degrees(151,12,30) j02 = [bb_1x, bb_1y] j12 = [bb_2x, bb_2y] j22 = [bb_3x, bb_3y] j32 = [bb_4x, bb_4y] j42 = [bb_5x, bb_5y] j52 = [bb_6x, bb_6y] j62 = [bb_7x, bb_7y] j72 = [bb_8x, bb_8y] j82 = [bb_9x, bb_9y] j92 = [bb_10x, bb_10y] botanybay_polygon_2, zone = convert_from_latlon_to_utm([j92, j12, j22, j62, j82, j72, j42]) # worked # close to botany bay opening (340000,6236000) # x0 = 25964 # y0 = 11049 # around 10km from botany bay opening (350000,6236000) # x0 = 35964 # y0 = 11049 # around 21km from botany bay opening (361000,6236000) #x0 = 46964 #y0 = 11049 # not used for sydney scenario, original interior regions listed though # setting up problem area for doing just around the harbour hsouth = degminsec2decimal_degrees(-33,54,0) hnorth = degminsec2decimal_degrees(-33,48,0) hwest = degminsec2decimal_degrees(151,0,0) heast = degminsec2decimal_degrees(151,30,0) hp0 = [hsouth, hwest] hp1 = [hsouth, heast] hp2 = [hnorth, heast] hp3 = [hnorth, hwest] polygon_h, zone = convert_from_latlon_to_utm([hp0, hp1, hp2, hp3]) #Interior regions - the Harbour - take 1 harbour_south = degminsec2decimal_degrees(-33,53,0) harbour_north = degminsec2decimal_degrees(-33,47,0) harbour_west = degminsec2decimal_degrees(151,5,0) harbour_east = degminsec2decimal_degrees(151,19,0) #harbour_south1 = degminsec2decimal_degrees(-33,53,0) #harbour_south2 = degminsec2decimal_degrees(-33,52,0) #harbour_north1 = degminsec2decimal_degrees(-33,45,0) #harbour_north2 = degminsec2decimal_degrees(-33,48,0) #harbour_west = degminsec2decimal_degrees(151,5,0) #harbour_east = degminsec2decimal_degrees(151,19,0) k0 = [harbour_south, harbour_west] k1 = [harbour_south, harbour_east] k2 = [harbour_north, harbour_east] k3 = [harbour_north, harbour_west] harbour_polygon, zone = convert_from_latlon_to_utm([k0, k1, k2, k3]) # setting up problem area for doing just around Botany Bay bsouth = degminsec2decimal_degrees(-33,56,0) bnorth = degminsec2decimal_degrees(-34,3,0) bwest = degminsec2decimal_degrees(151,0,0) beast = degminsec2decimal_degrees(151,30,0) bp0 = [bsouth, bwest] bp1 = [bsouth, beast] bp2 = [bnorth, beast] bp3 = [bnorth, bwest] polygon_bb, zone = convert_from_latlon_to_utm([bp0, bp1, bp2, bp3]) #Interior region - Botany Bay - take 1 botanybay_south = degminsec2decimal_degrees(-33,58,0) botanybay_north = degminsec2decimal_degrees(-34,1,0) botanybay_west = degminsec2decimal_degrees(151,5,0) botanybay_east = degminsec2decimal_degrees(151,18,0) j0 = [botanybay_south, botanybay_west] j1 = [botanybay_south, botanybay_east] j2 = [botanybay_north, botanybay_east] j3 = [botanybay_north, botanybay_west] botanybay_polygon, zone = convert_from_latlon_to_utm([j0, j1, j2, j3]) assert zone == refzone #x0 = 28964 + 42000 #y0 = 30049 #slump_origin = [x0+314036.58727982, y0+6224951.2960092] #Absolute UTM slump_origin = [385000.0, 6255000.0] #Absolute UTM # move 10km west slump_origin2 = [375000.0, 6255000.0] #Absolute UTM a = [340000, 6255000] b = [340000, 6270000] c = [318000, 6255000] d = [318000, 6270000] e = [316000, 6255000] f = [395000, 6255000] g = [355000, 6280000] h = [355000, 6224000] test_pts = [a, b, c, d, e, f, g, h] test_elev = [1.0, 4.0, 3.0, 0.1, 5, -10.0, -20, -15]