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

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

from os import sep, environ, getenv, getcwd
from os.path import expanduser
import sys
from time import localtime, strftime, gmtime
from anuga.utilities.polygon import read_polygon, plot_polygons, polygon_area, is_inside_polygon, number_mesh_triangles
from anuga.coordinate_transforms.redfearn import degminsec2decimal_degrees, convert_points_from_latlon_to_utm
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'

#time stuff
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
build_time = time+'_build'
run_time = time+'_run'
print 'gtime: ', gtime

tide = 0.6

#Making assumptions about the location of scenario data
state = 'western_australia'
scenario_name = 'dampier'
scenario = 'dampier_tsunami_scenario_2006'

# onshore data provided by WA DLI
onshore_name = 'dampier_dli_ext' # original
#island
island_name = 'rott_dli_ext' # original
island_name1 = 'gard_dli_ext'
island_name2 = 'carnac_island_dted'
island_name3 = 'penguin_dted'

# AHO + DPI data + colin French coastline
coast_name = 'waterline'
offshore_name = 'dampier_bathymetry'
offshore1_name = 'missing_fairsheets'

#final topo name
combined_name ='dampier_combined_elevation'
combined_smaller_name = 'dampier_combined_elevation_smaller'

topographies_in_dir = home+sep+state+sep+scenario+sep+'elevation_final'+sep+'points'+sep
topographies_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'topographies'+sep
topographies_time_dir = topographies_dir+build_time+sep

# input topo file location
onshore_in_dir_name = topographies_in_dir + onshore_name
island_in_dir_name = topographies_in_dir + island_name
island_in_dir_name1 = topographies_in_dir + island_name1
island_in_dir_name2 = topographies_in_dir + island_name2
island_in_dir_name3 = topographies_in_dir + island_name3

coast_in_dir_name = topographies_in_dir + coast_name
offshore_in_dir_name = topographies_in_dir + offshore_name
offshore1_in_dir_name = topographies_in_dir + offshore1_name

onshore_dir_name = topographies_dir + onshore_name
island_dir_name = topographies_dir + island_name
island_dir_name1 = topographies_dir + island_name1
island_dir_name2 = topographies_dir + island_name2
island_dir_name3 = topographies_dir + island_name3

coast_dir_name = topographies_dir + coast_name
offshore_dir_name = topographies_dir + offshore_name

#final topo files
combined_dir_name = topographies_dir + combined_name
combined_time_dir_name = topographies_time_dir + combined_name
combined_smaller_name_dir = topographies_dir + combined_smaller_name
#combined_time_dir_final_name = topographies_time_dir + combined_final_name

meshes_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'meshes'+sep
meshes_dir_name = meshes_dir + scenario_name

polygons_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'polygons'+sep
tide_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'tide_data'+sep

boundaries_source = '????'
#boundaries locations
boundaries_in_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'boundaries'+sep+'urs'+sep+boundaries_source+sep
boundaries_in_dir_name = boundaries_in_dir + scenario_name
boundaries_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'boundaries'+sep
boundaries_dir_name = boundaries_dir + scenario_name
#boundaries_time_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'boundaries'+sep+build_time+sep
#boundaries_time_dir_name = boundaries_time_dir + boundaries_name  #Used by post processing

#output locations
output_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'outputs'+sep
output_build_time_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'outputs'+sep+build_time+sep
output_run_time_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'outputs'+sep+run_time+sep
output_run_time_dir_name = output_run_time_dir + scenario_name  #Used by post processing

#gauges
gauge_name = 'dampier.csv'
gauges_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'gauges'+sep
gauges_dir_name = gauges_dir + gauge_name


###############################
# Domain definitions
###############################

refzone = 50 
south = degminsec2decimal_degrees(-20,55,0)
north = degminsec2decimal_degrees(-20,15,0)
west = degminsec2decimal_degrees(116,17,0)
east = degminsec2decimal_degrees(117,10,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]

poly_all, zone = convert_from_latlon_to_utm([p0, p1, p2, p3, p4, p5, p6, p7, p8])
refzone = zone
print 'Area of bounding polygon', polygon_area(poly_all)/1000000.0

res_poly_all = 100000

###############################
# Interior region definitions
###############################

poly_pos20_neg20 = read_polygon(polygons_dir+'pos20_neg20_pts.csv')
res_pos20_neg20 = 20000

poly_dampier = read_polygon(polygons_dir+'dampier_pts.csv')
res_dampier = 500

poly_karratha = read_polygon(polygons_dir+'karratha_pts.csv')
res_karratha = 500

poly_delambre = read_polygon(polygons_dir+'delambre_pts.csv')
res_delambre = 1000

poly_mainisland = read_polygon(polygons_dir+'mainisland_pts.csv')
res_mainisland = 1000

poly_NWislands = read_polygon(polygons_dir+'NWislands_pts.csv')
res_NWislands = 1000

plot_polygons([poly_pos20_neg20,poly_dampier,poly_karratha,poly_delambre,polylmainisland,
               polyNWislands,poly_all],output_run_time_dir + 'poly_pic')

interior_regions = [[poly_pos20_neg20,res_pos20_neg20],[poly_dampier,res_dampier],
                    [poly_karratha,res_karratha],[poly_delambre,res_delambre],
                    [poly_mainisland,res_mainisland],[poly_NWislands,res_NWislands]]
                    
trigs_min = number_mesh_triangles(interior_regions, poly_all, res_poly_all)

print 'min number triangles', trigs_min

###################################################################
# Clipping regions for export to asc and regions for clipping data
###################################################################

poly_bathy = read_polygon(polygons_dir+'polybathy.csv')

# exporting asc grid - Dampier
e_min_area = 474000
e_max_area = 480000
n_min_area = 7719000
n_max_area = 7725000

# exporting asc grid - Karratha
e_min_area = 
e_max_area = 
n_min_area = 
n_max_area = 

"""
# used in the CIPMA 2006 scenario
# 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

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')

# exporting asc grid
e_min_area = 474000
e_max_area = 480000
n_min_area = 7719000
n_max_area = 7725000

# used in the original 2005 scenario
#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
"""