source: anuga_work/production/dampier_2006/project.py @ 4063

"""Common filenames and locations for topographic data, meshes and outputs.
"""

import sys
from os import sep, environ, getenv, getcwd
from os.path import expanduser, basename

from anuga.coordinate_transforms.redfearn import\
     degminsec2decimal_degrees,\
     convert_from_latlon_to_utm

from time import localtime, strftime, gmtime, ctime
from anuga.geospatial_data.geospatial_data import *
from anuga.utilities.polygon import read_polygon, plot_polygons, polygon_area
from anuga.utilities.system_tools import get_user_name

# file and system info
#---------------------------------
codename = 'project.py'

home = getenv('INUNDATIONHOME') #Sandpit's parent dir   
user = get_user_name()

# INUNDATIONHOME is the inundation directory, not the data directory.
home += sep +'data'
#----------------------------------
# Location and naming of scenario data
#----------------------------------
state = 'western_australia'
scenario_name = 'dampier_tsunami'
scenario_datas_name = 'dampier_tsunami_scenario_2006'  #name of the directory where the data is stored
#scenario_datas_name = 'karratha_tsunami_scenario_2005' # Tmp location

#mesh_name = 'elevation50m'
boundaries_name = 'dampier'
boundaries_source = 'mag_9_corrected'
#boundaries_source = 'test'

tide = 2.4
#tide = 0.0

# topography file names
onshore_name = 'dli_no_islands'
coast_name = 'DTED_05_Contour'
islands_name = 'dted_islands'
offshore_name = 'XY100003902'
offshore_name1 = 'XY100003903'
offshore_name2 = 'XY100003951'
offshore_name3 = 'XY100006321'
offshore_name4 = 'XY100011756'
offshore_name5 = 'XY100014243'
offshore_name6 = 'XY100014244'
offshore_name7 = 'XY100021081'
offshore_name8 = 'XY100021082'
offshore_name9 = 'XY100021083'
offshore_name10 = 'XY100021085'
offshore_name11 = 'XY100021086'
offshore_name12 = 'XY100026309'
offshore_name13 = 'XY100026338'
offshore_name14 = 'XYDM83'

offshore_old = 'elevation50m'

combined_name ='dampier_combined_elevation'
combined_final_name ='dampier_combined_elevation_final'

gauge_name = 'dampier_gauges_up2.csv'

#Derive subdirectories and filenames

meshes_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'meshes'+sep
topographies_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'topographies'+sep
gauges_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'gauges'+sep
polygons_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'polygons'+sep
boundaries_in_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'boundaries'+sep+'urs'+sep+boundaries_source+sep
#outputdir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'output'+sep
tide_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'tide_data'+sep

time = strftime('%Y%m%d_%H%M%S',localtime()) #gets time for new dir
gtime = strftime('%Y%m%d_%H%M%S',gmtime()) #gets time for new dir
#cctime = strftime('%Y%m%d_%H%M%S',ctime()) #gets time for new dir
build_time = time+'_build'
run_time = time+'_run'

print 'gtime: ', gtime
output_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'outputs'+sep
output_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'outputs'+sep
output_build_time_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'outputs'+sep+build_time+sep
output_run_time_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'outputs'+sep+run_time+sep
topographies_time_dir = topographies_dir+build_time+sep
boundaries_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'boundaries'+sep
boundaries_time_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'boundaries'+sep+build_time+sep
meshes_time_dir = meshes_dir+build_time+sep

#ideas
#boundaries_time_dir = boundaries_in_dir+'urs'+sep+boundaries_source+sep

gauge_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'gauges'+sep
gauge_filename = gauge_dir + 'dampier_gauges_up2.csv'

gauges_dir_name = gauges_dir + gauge_name

onshore_dir_name = topographies_dir + onshore_name
coast_dir_name = topographies_dir + coast_name
islands_dir_name = topographies_dir + islands_name
offshore_dir_name = topographies_dir + offshore_name
offshore_dir_name1 = topographies_dir + offshore_name1
offshore_dir_name2 = topographies_dir + offshore_name2
offshore_dir_name3 = topographies_dir + offshore_name3
offshore_dir_name4 = topographies_dir + offshore_name4
offshore_dir_name5 = topographies_dir + offshore_name5
offshore_dir_name6 = topographies_dir + offshore_name6
offshore_dir_name7 = topographies_dir + offshore_name7
offshore_dir_name8 = topographies_dir + offshore_name8
offshore_dir_name9 = topographies_dir + offshore_name9
offshore_dir_name10 = topographies_dir + offshore_name10
offshore_dir_name11 = topographies_dir + offshore_name11
offshore_dir_name12 = topographies_dir + offshore_name12
offshore_dir_name13 = topographies_dir + offshore_name13
offshore_dir_name14 = topographies_dir + offshore_name14

offshore_dir_name_old = topographies_dir + offshore_old



#output dir
combined_dir_name = topographies_dir + combined_name
combined_time_dir_name = topographies_time_dir + combined_name
combined_time_dir_final_name = topographies_time_dir + combined_final_name


meshes_dir_name = meshes_dir + scenario_name
meshes_time_dir_name = meshes_time_dir + scenario_name
#output_build_time_dir_name = output_build_time_dir + scenario_name  #Used by post processing
output_run_time_dir_name = output_run_time_dir + scenario_name  #Used by post processing
boundaries_in_dir_name = boundaries_in_dir + boundaries_name
boundaries_time_dir_name = boundaries_time_dir + boundaries_name  #Used by post processing
boundaries_dir_name = boundaries_dir + boundaries_name



# Regions

refzone = 50 
south = degminsec2decimal_degrees(-20,55,0)
north = degminsec2decimal_degrees(-20,15,0)
#north = degminsec2decimal_degrees(-19,15,0)
west = degminsec2decimal_degrees(116,17,0)
east = degminsec2decimal_degrees(117,10,0)

#only used to clip boundary condition
#north_boundary = north + 0.02
#south_boundary = south - 0.02
#west_boundary = west - 0.02
#east_boundary = east + 0.02

south_boundary = degminsec2decimal_degrees(-21,0,0)
#north_boundary = degminsec2decimal_degrees(-19,00,0)
north_boundary = degminsec2decimal_degrees(-20,10,0)
#west_boundary = degminsec2decimal_degrees(116,0,0)
#east_boundary = degminsec2decimal_degrees(118,00,0)
west_boundary = degminsec2decimal_degrees(116,10,0)
east_boundary = degminsec2decimal_degrees(117,15,0)


p0 = [south, degminsec2decimal_degrees(116,32,0)]
p1 = [south, west]
p2 = [degminsec2decimal_degrees(-20,23,0), west]
p3 = [north, degminsec2decimal_degrees(116,45,0)]
p4 = [north, degminsec2decimal_degrees(117,0,0)]
p5 = [p2[0], degminsec2decimal_degrees(117,8,0)]
p6 = [degminsec2decimal_degrees(-20,30,0), east]
p7 = [degminsec2decimal_degrees(-20,38,0), east]
p8 = [south, east]

bounding_polygon, zone =\
                  convert_from_latlon_to_utm([p0, p1, p2, p3, p4, p5, p6, p7, p8])
#bounding_polygon, zone =\
#                  convert_from_latlon_to_utm([p1, p2, p3, p4, p5, p6, p7])
refzone = zone
from anuga.utilities.polygon import read_polygon, plot_polygons, polygon_area, is_inside_polygon
print 'poly area', polygon_area(bounding_polygon)/1000000.0

#Interior regions

# CIPMA point of interest
cipma_latitude = -20.588456
cipma_longitude = 116.771527


k0 = [cipma_latitude-0.02, cipma_longitude-0.02]
k1 = [cipma_latitude-0.02, cipma_longitude+0.02]
k2 = [cipma_latitude+0.02, cipma_longitude+0.02]
k3 = [cipma_latitude+0.02, cipma_longitude-0.02]

cipma_polygon, zone = convert_from_latlon_to_utm([k0, k1, k2, k3])
assert zone == refzone

e_min_area = 474000
e_max_area = 480000
n_min_area = 7719000
n_max_area = 7725000

e_min_area_mom = 442200
e_max_area_mom = 497028
n_min_area_mom = 7720095
n_max_area_mom = 7750500

poly_facility = read_polygon(polygons_dir+'facility.csv')

poly_pipeline = read_polygon(polygons_dir+'pipeline2.csv')

poly_interior = read_polygon(polygons_dir+'interior.csv')

poly_coast = read_polygon(polygons_dir+'coast_final.csv')

clip_poly_e = read_polygon(polygons_dir+'gap_e.csv')

clip_poly_nw = read_polygon(polygons_dir+'gap_nw.csv')

clip_poly_mid_w = read_polygon(polygons_dir+'gap_mid_w.cvs')

clip_poly_mid_e = read_polygon(polygons_dir+'gap_mid_e.cvs')

plot_polygons([bounding_polygon,poly_facility,poly_pipeline,poly_interior,poly_coast],'polys',verbose=True)
#Interior regions
karratha_south = degminsec2decimal_degrees(-20,44,0)
karratha_north = degminsec2decimal_degrees(-20,42,0)
karratha_west = degminsec2decimal_degrees(116,48,0)
karratha_east = degminsec2decimal_degrees(116,53,30)

k0 = [karratha_south, karratha_west]
k1 = [karratha_south, karratha_east]
k2 = [karratha_north, karratha_east]
k3 = [karratha_north, karratha_west]    

karratha_polygon, zone = convert_from_latlon_to_utm([k0, k1, k2, k3])
assert zone == refzone


#Interior regions
dampier_south = degminsec2decimal_degrees(-20,40,0)
dampier_north = degminsec2decimal_degrees(-20,38,10)
dampier_west = degminsec2decimal_degrees(116,43,0)
dampier_east = degminsec2decimal_degrees(116,45,0)

d0 = [dampier_south, dampier_west]
d1 = [dampier_south, dampier_east]
d2 = [dampier_north, dampier_east]
d3 = [dampier_north, dampier_west]    

dampier_polygon, zone = convert_from_latlon_to_utm([d0, d1, d2, d3])
assert zone == refzone


#Interior regions
refinery_south = degminsec2decimal_degrees(-20,37,50)
refinery_north = degminsec2decimal_degrees(-20,36,0)
refinery_west = degminsec2decimal_degrees(116,44,0)
refinery_east = degminsec2decimal_degrees(116,46,10)

d0 = [refinery_south, refinery_west]
d1 = [refinery_south, refinery_east]
d2 = [refinery_north, refinery_east]
d3 = [refinery_north, refinery_west]    

refinery_polygon, zone = convert_from_latlon_to_utm([d0, d1, d2, d3])
assert zone == refzone


#Interior region around 468899, 7715177:
#lat (-20, 39, 44.93753), lon (116, 42, 5.09106)

point_south = degminsec2decimal_degrees(-20,39,46)
point_north = degminsec2decimal_degrees(-20,39,42)
point_west = degminsec2decimal_degrees(116,42,0)
point_east = degminsec2decimal_degrees(116,42,10)

d0 = [point_south, point_west]
d1 = [point_south, point_east]
d2 = [point_north, point_east]
d3 = [point_north, point_west]    

point_polygon, zone = convert_from_latlon_to_utm([d0, d1, d2, d3])
assert zone == refzone


#Neils areas around interesting points
neil1_point1 = [degminsec2decimal_degrees(-20,35,34),
                degminsec2decimal_degrees(116,45,18)]
neil1_point2 = [degminsec2decimal_degrees(-20,36,15),
                degminsec2decimal_degrees(116,46,18)]
neil1_point3 = [degminsec2decimal_degrees(-20,35,9),
                degminsec2decimal_degrees(116,47,17)]
neil1_point4 = [degminsec2decimal_degrees(-20,34,26),
                degminsec2decimal_degrees(116,46,17)]

neil1_polygon, zone = convert_from_latlon_to_utm([neil1_point1,
                                                  neil1_point2,
                                                  neil1_point3,
                                                  neil1_point4])
assert zone == refzone




neil2_point1 = [degminsec2decimal_degrees(-20,39,36),
                degminsec2decimal_degrees(116,41,33)]
neil2_point2 = [degminsec2decimal_degrees(-20,40,10),
                degminsec2decimal_degrees(116,42,13)]
neil2_point3 = [degminsec2decimal_degrees(-20,38,39),
                degminsec2decimal_degrees(116,43,49)]
neil2_point4 = [degminsec2decimal_degrees(-20,38,5),
                degminsec2decimal_degrees(116,43,9)]

neil2_polygon, zone = convert_from_latlon_to_utm([neil2_point1,
                                                  neil2_point2,
                                                  neil2_point3,
                                                  neil2_point4])
assert zone == refzone





#Withnell bay
wb_point1 = [degminsec2decimal_degrees(-20,35,34),
                degminsec2decimal_degrees(116,45,18)]
wb_point2 = [degminsec2decimal_degrees(-20,36,15),
                degminsec2decimal_degrees(116,46,18)]
wb_point3 = [degminsec2decimal_degrees(-20,35,9),
                degminsec2decimal_degrees(116,47,17)]
wb_point4 = [degminsec2decimal_degrees(-20,34,26),
                degminsec2decimal_degrees(116,46,17)]

wb_polygon, zone = convert_from_latlon_to_utm([wb_point1, wb_point2,
                                               wb_point3, wb_point4])
assert zone == refzone





#Larger Withnell bay
lwb_point1 = [degminsec2decimal_degrees(-20,35,59),
                degminsec2decimal_degrees(116,42,00)]
lwb_point2 = [degminsec2decimal_degrees(-20,36,50),
                degminsec2decimal_degrees(116,46,50)]
lwb_point3 = [degminsec2decimal_degrees(-20,34,00),
                degminsec2decimal_degrees(116,47,39)]
lwb_point4 = [degminsec2decimal_degrees(-20,33,00),
                degminsec2decimal_degrees(116,42,50)]

lwb_polygon, zone = convert_from_latlon_to_utm([lwb_point1, lwb_point2,
                                                lwb_point3, lwb_point4])
                                                      
assert zone == refzone