source: anuga_work/production/sydney_2006/project.py @ 5001

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