source: production/sydney_2006/project.py @ 2763

"""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]