source: anuga_work/development/convergence_okushiri_2008/project_truescale.py @ 6409

"""Common filenames and run parameters for truescale Okushiri Island
convergence study. Formats are given as ANUGA native netCDF where applicable.

"""
from os import sep, environ, getenv, getcwd, umask
from os.path import expanduser, basename, join
from anuga.utilities.polygon import read_polygon, plot_polygons, polygon_area, is_inside_polygon, number_mesh_triangles
import sys
from anuga.coordinate_transforms.redfearn import degminsec2decimal_degrees
from time import localtime, strftime, gmtime
from anuga.utilities.system_tools import get_user_name, get_host_name

home = join(getenv('INUNDATIONHOME'),'data', 'anuga_validation',
            'convergence_okushiri_2008') # Location of Data
user = get_user_name()
host = get_host_name()
#needed when running using mpirun, mpirun doesn't inherit umask from .bashrc
umask(002)

#-------------------
# Input file names 
#-------------------

# Given boundary wave
boundary_filename = 'okushiri_truescale_input.tms'

# Observed timeseries
validation_filename = 'okushiri_output_truescale_ch5-7-9.txt'

# Digital Elevation Model
bathymetry_filename = 'okushiri_truescale_bathymetry.pts'


#------------------------------------
# Output file names and directories 
#------------------------------------

# Model output
output_filename = 'okushiri_truescale.sww'

# Time stuff
time = strftime('%Y%m%d_%H%M%S',gmtime()) #gets time for new dir
run_time = time+'_run'

# Set anuga input directory names

anuga_dir = join(home,'anuga')+sep

mesh_dir = join(anuga_dir, 'meshes')+sep
mesh_name = join(mesh_dir, 'okushiri_truescale')

polygons_dir = join(anuga_dir, 'polygons')+sep # Created with ArcGIS (csv files)

#------------------------
# Run parameters
#------------------------

finaltime=450
setup='sixteen'
polygons = 'contour_polygons'

if setup =='no polygons':
    print 'no interior polygons'
    base_resolution=1
    yieldstep=1
if setup =='octuple':
    print '8 times original resolution'
    base_resolution=0.125
    yieldstep=1
if setup =='sixteen':
    print '16 times original resolution'
    base_resolution=0.0625
    yieldstep=1
if setup =='quadruple':
    print '4 times original resolution'
    base_resolution=0.25
    yieldstep=1
if setup =='double':
    print 'double original resolution'
    base_resolution=0.5
    yieldstep=1
if setup =='1.5':
    print '1.5 times original resolution'
    base_resolution=0.66
    yieldstep=1
if setup =='original':
    print 'original resolution'
    base_resolution=1
    yieldstep=1
if setup =='0.75':
    print '0.75 times original resolution'
    base_resolution=1.33
    yieldstep=1
if setup =='half':
    print 'half original resolution'
    base_resolution=2
    yieldstep=1
if setup =='quarter':
    print '1/4 original resolution'
    base_resolution=4
    yieldstep=1
if setup =='eighth':
    print '1/8 original resolution'
    base_resolution=8
    yieldstep=1
if setup =='sixteenth':
    print '1/16 original resolution'
    base_resolution=16
    yieldstep = 1
if setup =='1-32':
    print '1/32 original resolution'
    base_resolution=32
    yieldstep=1
if setup =='1-64':
    print '1/64 original resolution'
    base_resolution=64
    yieldstep=1
  
  
#------------------------------
# Polygon definitions
#------------------------------
    
poly_all = read_polygon(polygons_dir+'bounding_polygon.csv')
res_poly_all = 16000*base_resolution

# Original polygon definitions

poly_gulleys = read_polygon(polygons_dir+'gulleys_polygon.csv')
res_gulleys = 3.2*base_resolution

poly_island = read_polygon(polygons_dir+'island_polygon.csv')
res_island = 32*base_resolution

poly_rhs = read_polygon(polygons_dir+'rhs_polygon.csv')
res_rhs = 80*base_resolution

                    
# Contour-based polygon definitions

poly_25 = read_polygon(polygons_dir+'polygon_25m.csv')
res_poly_25 = 800*base_resolution

poly_10 = read_polygon(polygons_dir+'polygon_10m.csv')
res_poly_10 = 80*base_resolution

poly_5 = read_polygon(polygons_dir+'polygon_5m.csv')
res_poly_5 = 32*base_resolution

poly_1 = read_polygon(polygons_dir+'polygon_1m.csv')
res_poly_1 = 10*base_resolution

                    
if polygons =='original_polygons':
    print 'original polygon definition'
    interior_regions = [[poly_gulleys,res_gulleys],[poly_island,res_island],
                    [poly_rhs,res_rhs]]

if polygons =='contour_polygons':
    print 'contour-based polygon definition'
    interior_regions = [[poly_25,res_poly_25],[poly_10,res_poly_10],
                    [poly_5,res_poly_5],[poly_1,res_poly_1] ]


trigs_min = number_mesh_triangles(interior_regions, poly_all, res_poly_all)
print 'min number triangles', trigs_min

this_area = polygon_area(poly_all)
print this_area


# Polygon QC
#polygon_plot='polygons.png'
#plot_polygons([poly_all, poly_25, poly_10, poly_5, poly_1], style=None, figname=polygon_plot)


#--------------------------------
# Set anuga output locations
#--------------------------------

# Directory comment
dir_comment='_'+setup+'_'+polygons+'_'+str(user)

output_dir = join(anuga_dir, 'outputs')+sep
output_run_time_dir = output_dir+run_time+dir_comment+sep
output_run_time_dir_name = output_run_time_dir + output_filename #Used by post processing

# Gauges
gauges_dir = join(anuga_dir,'gauges')+sep
gauge_name = 'gauge_location_okushiri.csv'
gauges_dir_name = gauges_dir+gauge_name

# Vertex coordinates
vertex_filename = output_run_time_dir+'vertex_coordinates.txt'


#--------------------------------
# Area definitions for sww2dem
#--------------------------------

xminDeep = 0
xmaxDeep = 950
yminDeep = 0
ymaxDeep = 1400

xminMid = 1000
xmaxMid = 1700
yminMid = 0
ymaxMid = 1400

xminShallow = 1725
xmaxShallow = 2175
yminShallow = 0
ymaxShallow = 1400