source: anuga_work/production/hobart_2006/project.py @ 3972

# -*- 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
from anuga.coordinate_transforms.redfearn import degminsec2decimal_degrees, convert_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 = 'tasmania'
scenario_dir_name = 'hobart_tsunami_scenario_2006'

# data provided by Tas SES and checked by NM&I
#onshore_name = 'hob3_topo' # original
onshore_name_25 = 'hob5_topo_25m' # 25m grid and clipped to 100m elevation or 3000m from coast
onshore_name = 'hob5_topo' # 12.5m grid and clipped to 100m elevation or 3000m from coast
#onshore_name_25 = 'hob6_topo_25m' # 25m grid NOT clipped
offshore_name_tas1 = 'derwent_2m'
offshore_name_tas2 = 'derwent_5m'
offshore_name_tas3 = 'south_east_tas' #actually this is AHO
offshore_name_tas4 = 'hobart_1m'

# AHO data and checked by NM&I
offshore_name1 = 'xy100003760'
offshore_name2 = 'xy100003761'
offshore_name3 = 'xy100003762'
offshore_name4 = 'xy100003907'
offshore_name5 = 'xy100003908'
offshore_name6 = 'xy100003909'
offshore_name7 = 'xy100003910'
offshore_name8 = 'xy100003932'
offshore_name9 = 'xy100003933'
offshore_name10 = 'xy100003934'
offshore_name11 = 'xy100003935'
offshore_name12 = 'xy100003936'
offshore_name13 = 'xy100003964'
offshore_name14 = 'xy100014250'
offshore_name15 = 'xy100014253'
offshore_name16 = 'xy100016142'

# developed by NM&I
coast_name = 'coastline_points'

boundary_basename = 'puysegur' # Mw 8.7
#boundary_basename = 'puysegur_clip' # Mw 8.5

#swollen/ all data output
basename = 'source'

codename = 'project.py'

#Derive subdirectories and filenames
#time = strftime('%Y%m%d_%H%M%S',localtime()) #gets time for new dir
local_time = strftime('%Y%m%d_%H%M%S',gmtime()) #gets time for new dir
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 + 'hobart_gauges_final.csv'
#gauge_filename = gaugedir + 'hobartGauges20061019.csv'
gauge_filename = gaugedir + 'gauges_MRT.csv'
gauge_filename_bom = gaugedir + 'gauges_bruny_fine.csv'
#gauge_filename = gaugedir + 'check_boundary.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
onshore_dem_name_25 = datadir + onshore_name_25
all_onshore_dem_name = datadir + 'combined_onshore'
offshore_dem_name_local1 = datadir + offshore_name_tas1
offshore_dem_name_local2 = datadir + offshore_name_tas2
offshore_dem_name_local3 = datadir + offshore_name_tas3
offshore_dem_name_local4 = datadir + offshore_name_tas4
offshore_dem_name_aho1 = datadir + offshore_name1
offshore_dem_name_aho2 = datadir + offshore_name2
offshore_dem_name_aho3 = datadir + offshore_name3
offshore_dem_name_aho4 = datadir + offshore_name4
offshore_dem_name_aho5 = datadir + offshore_name5
offshore_dem_name_aho6 = datadir + offshore_name6
offshore_dem_name_aho7 = datadir + offshore_name7
offshore_dem_name_aho8 = datadir + offshore_name8
offshore_dem_name_aho9 = datadir + offshore_name9
offshore_dem_name_aho10 = datadir + offshore_name10
offshore_dem_name_aho11 = datadir + offshore_name11
offshore_dem_name_aho12 = datadir + offshore_name12
offshore_dem_name_aho13 = datadir + offshore_name13
offshore_dem_name_aho14 = datadir + offshore_name14
offshore_dem_name_aho15 = datadir + offshore_name15
offshore_dem_name_aho16 = datadir + offshore_name16
coast_dem_name = datadir + coast_name

# addition once total grid delivered
onshore_offshore_dem_name_25 = datadir + '25m_se_tas' #25m grid
onshore_offshore_dem_name    = datadir + '50m_se_tas' #50m grid

# output names
bruny_dem_name_25   = datadir + 'bruny_25_dem'
hobart_dem_name_25  = datadir + 'hobart_25_dem'
combined_dem_name   = datadir + 'hobart_combined_elevation'
combined_dem_name_2 = datadir + 'hobart_combined_elevation_2'
combined_dem_name_3 = datadir + 'hobart_combined_elevation_jane'
combined_dem_name_4 = datadir + 'hobart_combined_elevation_test'

#outputname = outputtimedir + basename  #Used by post processing

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

# bounding box for clipping MOST output (much bigger than study area)
south = degminsec2decimal_degrees(-44,45,0)
north = degminsec2decimal_degrees(-42,0,0)
west  = degminsec2decimal_degrees(146,45,0)
east  = degminsec2decimal_degrees(148,25,0) 

###Main Domain of Hobart: 
d0 = [south, west]
d1 = [south, east]
d2 = [north, east]
d3 = [north, west]
polyAll2, zone = convert_from_latlon_to_utm([d0, d1, d2, d3])
refzone = zone

# Mark run - morning Fri 29 Sep; surrounds -100m and 20mish elevation
polyAll = read_polygon(polygondir+'new_extent_2.csv')
#plot_polygons([polyAll, polyAll2],'boundingpoly',verbose=False)
#print 'Area of bounding polygon', polygon_area(polyAll)
polyAll_refine = read_polygon(polygondir+'extent_refine.csv')
#print 'Area of refined bounding polygon', polygon_area(polyAll_refine)
#print 'change', polygon_area(polyAll_refine)/polygon_area(polyAll)

polyAll_refine2 = read_polygon(polygondir+'export_refine_2.csv')
#print 'Area of refined bounding polygon 2', polygon_area(polyAll_refine2)
#print 'change 2', polygon_area(polyAll_refine2)/polygon_area(polyAll)

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

# clipping 12.5m onshore data set
eastingmin = 520000
eastingmax = 536000
northingmin = 5245000
northingmax = 5260000
# for original model setup
# clipping 25m data set - Hobart
eastingmin25 = 524208.387
eastingmax25 = 554867.24
northingmin25 = 5229154.555
northingmax25 = 5258511.857
hob_clip = [[eastingmin25,northingmin25],[eastingmin25,northingmax25],[eastingmax25,northingmax25],[eastingmax25,northingmin25]]
# clipping 25m data set - Bruny
eastingmin25_2 = 523672.502
eastingmax25_2 = 536020.057
northingmin25_2 = 5203120.035
northingmax25_2 = 5212052.309
bruny_clip = [[eastingmin25_2,northingmin25_2],[eastingmin25_2,northingmax25_2],[eastingmax25_2,northingmax25_2],[eastingmax25_2,northingmin25_2]]
# overall 25m clipping region - covering Hobart and Bruny
eastingmin25_3 = 523672.502
eastingmax25_3 = 554867.24
northingmin25_3 = 5203120.035
northingmax25_3 = 5258511.857

# export for refined grids at paleo sites
eastingmin_site13 = 532835.94
eastingmax_site13 = 533848.23
northingmin_site13 = 5237384.18
northingmax_site13 = 5238195.03

eastingmin_kingston = 525500 #525905.49 
eastingmax_kingston = 527773.42
northingmin_kingston = 5240256.71
northingmax_kingston = 5242300 #5242148.58

eastingmin_bruny = 523638.9
eastingmax_bruny = 533000.88
northingmin_bruny = 5204967.41
northingmax_bruny = 5212177.69

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

# Hobart digitized polygons
poly_hobart1 = read_polygon(polygondir+'Hob_poly1.csv')
poly_hobart2 = read_polygon(polygondir+'Hob_poly2.csv')
poly_hobart3 = read_polygon(polygondir+'Hob_poly3.csv')
poly_hobart4 = read_polygon(polygondir+'Hob_poly4.csv')
poly_hobart5 = read_polygon(polygondir+'new_interior_region.csv')

#plot_polygons([polyAll, poly_hobart1,poly_hobart2,poly_hobart3,poly_hobart5],'orig_model',verbose=False)

# to refine around paleo sites
poly_kingston = read_polygon(polygondir+'Kingston.csv')
poly_bruny = read_polygon(polygondir+'bruny_refine.csv')
poly_site13 = read_polygon(polygondir+'Site13.csv')
#print 'areas of refined polys', polygon_area(poly_kingston), polygon_area(poly_bruny), polygon_area(poly_site13)

#plot_polygons([polyAll, poly_kingston,poly_bruny,poly_site13,poly_hobart5],'orig_model',verbose=False)
#print 'Area of interior regions', (polygon_area(poly_hobart1)+\
#      polygon_area(poly_hobart2)+polygon_area(poly_hobart3)+\
#      polygon_area(poly_hobart4)+polygon_area(poly_hobart5))/1000000.

#print 'Area of original 25m data clipping', (polygon_area(hob_clip)+polygon_area(bruny_clip))/1000000.
#print 'Area of new interior region', polygon_area(poly_hobart5)/1000000.