source: anuga_work/production/broome_2006/project.py @ 3788

# -*- coding: cp1252 -*-
"""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
from anuga.coordinate_transforms.redfearn import degminsec2decimal_degrees, convert_points_from_latlon_to_utm

if sys.platform == 'win32':
    home = getenv('INUNDATIONHOME')
    user = getenv('USERPROFILE')

else:    
    home = getenv('INUNDATIONHOME', sep+'d'+sep+'xrd'+sep+'gem'+sep+'2'+sep+'ramp'+sep+'risk_assessment_methods_project'+sep+'inundation')     
    user = getenv('LOGNAME')
    print 'USER:', user

# INUNDATIONHOME is the inundation directory, not the data directory.
home += sep +'data'

#Making assumptions about the location of scenario data
state = 'western_australia'
scenario_dir_name = 'broome_tsunami_scenario_2006'

# onshore data provided by WA DLI
onshore_name = '' # original

# offshore data provided by WA DPI
offshore_name_dpi1 = ''

# AHO data
offshore_name1 = 'xy100003760'

# developed by NM&I
coast_name = 'coastline_points'

boundary_basename = 'SU-AU' # Mw ?

#swollen/ all data output
basename = 'source'
codename = 'project.py'

#Derive subdirectories and filenames
local_time = strftime('%Y%m%d_%H%M%S',gmtime()) 
meshdir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'meshes'+sep
datadir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'topographies'+sep
gaugedir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'gauges'+sep
polygondir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'polygons'+sep
boundarydir = home+sep+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'boundaries'+sep
outputdir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'outputs'+sep
outputtimedir = outputdir + local_time + sep
polygondir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'polygons'+sep

gauge_filename = gaugedir + 'broome_gauges.csv'

codedir = getcwd()+sep                            
codedirname = codedir + 'project.py'
meshname = outputtimedir + 'mesh_' + basename

# Necessary if using point datasets, rather than grid
onshore_dem_name = datadir + onshore_name
offshore_dem_name_local1 = datadir + offshore_name_dpi1
offshore_dem_name_aho1 = datadir + offshore_name1

coast_dem_name = datadir + coast_name

combined_dem_name   = datadir + 'broome_combined_elevation'

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

# bounding box for clipping MOST/URS output (much bigger than study area)
south = degminsec2decimal_degrees(-19,0,0)
north = degminsec2decimal_degrees(-17,15,0)
west  = degminsec2decimal_degrees(120,0,0)
east  = degminsec2decimal_degrees(122,30,0) 

d0 = [south, west]
d1 = [south, east]
d2 = [north, east]
d3 = [north, west]
poly_bc, zone = convert_points_from_latlon_to_utm([d0, d1, d2, d3])
refzone = zone

# bounding polygon for study area
polyAll = read_polygon(polygondir+'extent.csv')

# plot bounding polygon and make sure BC info surrounds it
plot_polygons([polyAll, poly_bc],'boundingpoly',verbose=False)
print 'Area of bounding polygon', polygon_area(polyAll)/1000000.0

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

# clipping 20m onshore data set
#eastingmin = 520000
#eastingmax = 536000
#northingmin = 5245000
#northingmax = 5260000

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

# broome digitized polygons
poly_broome1 = read_polygon(polygondir+'Broome_Local_Polygon.csv')
poly_broome2 = read_polygon(polygondir+'Broome_Close2.csv')
poly_broome3 = read_polygon(polygondir+'Broome_Coast.csv')
poly_broome4 = read_polygon(polygondir+'Cable_Beach.csv')

plot_polygons([polyAll,poly_broome1,poly_broome2,poly_broome3,poly_broome4],'boundingpoly2',verbose=False)
print 'Area of local polygon', polygon_area(poly_broome1)/1000000.0
print 'Area of close polygon', polygon_area(poly_broome2)/1000000.0
print 'Area of coastal polygon', polygon_area(poly_broome3)/1000000.0
print 'Area of cable beach polygon', polygon_area(poly_broome4)/1000000.0

for i in poly_broome3:
    v = is_inside_polygon(i,poly_broome1, verbose=False)
    if v == False: print v