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

# -*- 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 = 'dted2_z51' # original

# AHO + DPI data
offshore_name1 = 'XY100011610'
offshore_name2 = 'XY100011611'
offshore_name3 = 'XY100011613'
offshore_name4 = 'XY100011614'
offshore_name5 = 'XY100011616'
offshore_name6 = 'XY100011617'
offshore_name7 = 'XY100011618'
offshore_name8 = 'XY100011621'
offshore_name9 = 'XY100011623'
offshore_name10 = 'XY100011745'
offshore_name11 = 'XY100011746'
offshore_name12 = 'XY100017530'
offshore_name13 = 'XY100017532'
offshore_name14 = 'XY100017538'
offshore_name15 = 'XY100017540'
offshore_name16 = 'XYBR66'
offshore_name17 = 'XYBR70'
offshore_name18 = 'XYBR80'
offshore_name19 = 'XYBR88'
offshore_name20 = 'XYBR93'
offshore_name21 = 'XYBR0110'
offshore_name22 = 'XYWADPI'

# developed by NM&I
coast_name = 'coastline'

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'
community_filename = gaugedir + 'CHINS_v2.csv'
community_broome = gaugedir + 'community_broome.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_name1 = datadir + offshore_name1
offshore_dem_name2 = datadir + offshore_name2
offshore_dem_name3 = datadir + offshore_name3
offshore_dem_name4 = datadir + offshore_name4
offshore_dem_name5 = datadir + offshore_name5
offshore_dem_name6 = datadir + offshore_name6
offshore_dem_name7 = datadir + offshore_name7
offshore_dem_name8 = datadir + offshore_name8
offshore_dem_name9 = datadir + offshore_name9
offshore_dem_name10 = datadir + offshore_name10
offshore_dem_name11 = datadir + offshore_name11
offshore_dem_name12 = datadir + offshore_name12
offshore_dem_name13 = datadir + offshore_name13
offshore_dem_name14 = datadir + offshore_name14
offshore_dem_name15 = datadir + offshore_name15
offshore_dem_name16 = datadir + offshore_name16
offshore_dem_name17 = datadir + offshore_name17
offshore_dem_name18 = datadir + offshore_name18
offshore_dem_name19 = datadir + offshore_name19
offshore_dem_name20 = datadir + offshore_name20
offshore_dem_name21 = datadir + offshore_name21
offshore_dem_name22 = datadir + offshore_name22

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
###################################################################

# exporting asc grid
eastingmin = 406215.87
eastingmax = 440208.78
northingmin = 7983427.73
northingmax = 8032834.52

        
        


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

# broome digitized polygons
poly_broome1 = read_polygon(polygondir+'Broome_Local_Polygon_update.csv')
poly_broome2 = read_polygon(polygondir+'Broome_Close2_update.csv')
poly_broome3 = read_polygon(polygondir+'Broome_Coast_update.csv')
#poly_broome4 = read_polygon(polygondir+'Cable_Beach_revised.csv')

plot_polygons([polyAll,poly_broome1,poly_broome2,poly_broome3],'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

def number_mesh_triangles(interior_regions, bounding_poly, remainder_res):
    from anuga.utilities.polygon import polygon_area
    
    # TO DO check if any of the regions fall inside one another
    no_triangles = 0.0
    area = polygon_area(bounding_poly)
    for i,j in interior_regions:
        this_area = polygon_area(i)
        no_triangles += this_area/j
        area -= this_area
        print j, this_area/1000000., area/1000000.
    no_triangles += area/remainder_res
    return int(no_triangles/0.7)