"""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 utilities.polygon import read_polygon import sys from pmesh.create_mesh import convert_points_from_latlon_to_utm #Making assumptions about the location of scenario data 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_friction04' basename4 = 'slump_poly_ingo_test' if sys.platform == 'win32': home = environ['INUNDATIONHOME'] #Sandpit's parent dir else: home = expanduser('~') #Derive subdirectories and filenames meshdir = home+sep+scenario_dir_name+sep+'meshes'+sep datadir = home+sep+scenario_dir_name+sep+'topographies'+sep outputdir = home+sep+scenario_dir_name+sep+'output'+sep polygondir = home+sep+scenario_dir_name+sep+'polygons'+sep gaugedir = home+sep+scenario_dir_name+sep+'gauges'+sep meshname = meshdir + basename meshname4 = meshdir + basename4 coarsedemname = datadir + coarsename finedemname = datadir + finename combineddemname = datadir + 'sydneytopo' outputname = outputdir + basename #Used by post processing outputname4 = outputdir + basename4 #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' 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_Benfield' polygonptsfile = polygondir + 'poly' integraltimeseries = outputdir + 'integral_time_series' #Georeferencing from 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_points_from_latlon_to_utm([p0, p1, p2, p3]) refzone = zone print 'Got refzone', refzone dsouth = degminsec2decimal_degrees(-34,05,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,2,20) 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] diffpolygonall, zone = convert_points_from_latlon_to_utm([dp0, dp1, dp2, dp3, dp4, dp5, dp6, dp7]) # 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] # didn't like this one - poly2 from Ingo # warning generated about vertex blah doesn't belong to an element # is there a maximum number of vertices per polygon? p0 = [343488.9455, 6270644.956] p1 = [343681.208, 6270500.759] p2 = [343997.0678, 6270253.564] p3 = [343482.0789, 6269930.838] p4 = [343550.7441, 6269546.313] p5 = [343111.287, 6269099.989] p6 = [342630.6307, 6267534.423] p7 = [342870.9588, 6266380.848] p8 = [342479.5673, 6266126.787] p9 = [342174.4759, 6266317.303] p10 = [342248.668, 6266727.125] p11 = [342101.9088, 6267108.699] p12 = [340986.2201, 6267454.634] p13 = [340663.6936, 6266113.602] p14 = [340120.4909, 6265892.926] p15 = [339560.3132, 6266317.303] p16 = [339237.7867, 6266928.406] p17 = [339000.1355, 6267912.961] p18 = [339169.8863, 6268642.889] p19 = [340001.6654, 6269016.341] p20 = [340392.0922, 6269423.743] p21 = [340867.3945, 6268863.566] p22 = [341614.2982, 6268659.865] p23 = [341495.4726, 6269220.042] p24 = [341750.0988, 6269287.943] p25 = [341902.8746, 6269882.07] p26 = [342700.7034, 6269780.22] p27 = [343125.0805, 6269746.27] p28 = [343125.0805, 6270221.572] p29 = [343057.1801, 6270527.124] p30 = [343488.9455, 6270644.956] testpoly = [p0,p1,p2,p3,p4,p5,p6,p7,p8,p9,p10,p11,p12,p13,p14,p15,p16,p17,p18,p19,p20,p21,p22,p23,p24,p25,p26,p27,p28,p29,p30] # 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] 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] #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] 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] 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] #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_points_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_points_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_points_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_points_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_points_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_points_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 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]