source: anuga_work/production/shark_bay_2007/project.py @ 4856

"""Common filenames and locations for topographic data, meshes and outputs.
Also includes origin for slump scenario.

All filenames are given without extension
"""

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

codename = 'project.py' # FIXME can be obtained automatically as __file__

home = join(getenv('INUNDATIONHOME'), 'data') # 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)

#Making assumptions about the location of scenario data
state = 'western_australia'
scenario_name = 'shark_bay'
scenario = 'shark_bay_tsunami_scenario'

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

#tide = 0.0 # MSL
tide = 0.79 # HAT average between Denham and Canarvon
#tide = -0.80 # Estimated LAT

# For frequency response study
amplitude = 0.5
period = 3000

#momentum_scale=50 # Experiment
momentum_scale=1 # Experiment                        

#Maybe will try to make project a class to allow these parameters to be passed in.
alpha = 0.1
friction=0.01
finaltime=40000
starttime=8000 # Action starts around 9000
setup='final'
#setup='trial'
source='shark_bay'

boundary_event = 'july2006'
#boundary_event = '10000'
#boundary_event = 'experimental'

if setup =='trial':
    print'trial'
    res_factor=10
    time_thinning=48
    yieldstep=240
if setup =='basic': 
    print'basic'
    res_factor=4
    time_thinning=12
    yieldstep=120
if setup =='final': 
    print'final'
    res_factor=1.0
    time_thinning=1
    yieldstep=10



#dir_comment='_'+setup+'_'+str(tide)+'_'+str(source)+'_'+\
#             str(user)+'_'+boundary_event+'_X'+str(momentum_scale)+'_fixedbathy'

dir_comment='_'+setup+'_'+str(tide)+'_'+str(source)+'_'+\
             str(user)+'_'+boundary_event


if boundary_event == 'experimental':
    dir_comment += '_A'+str(amplitude)+'_T'+str(period)
    

# elevation data filenames
ascii_grid_filenames = ['10m_dem_without_survey', '50m_dem_without_10m_dem',
                        'bathysteeppt', 'bathyleft', 'bathyright']
point_filenames = ['field_survey_north', 'field_survey_south',
                   'clipped_bathymetry_final', 'coast_points_final']

# final topo name
combined_name ='shark_bay_combined_elevation'
combined_small_name = 'shark_bay_combined_elevation_small'


anuga_dir = join(home, state, scenario, 'anuga')

topographies_in_dir = join(home, state, scenario,
                           'elevation_final', 'test_area')

topographies_dir = join(anuga_dir, 'topographies') # Output dir for ANUGA

# Convert topo file locations into absolute pathnames
for filename_list in [ascii_grid_filenames, point_filenames]:
    for i, filename in enumerate(filename_list):
        filename_list[i] = join(topographies_in_dir, filename)

#final topo files
combined_dir_name = join(topographies_dir, combined_name)
combined_small_dir_name = join(topographies_dir, combined_small_name)


mesh_dir = join(anuga_dir, 'meshes')
mesh_name = join(mesh_dir, scenario_name)


polygons_dir = join(anuga_dir, 'polygons') # Created with ArcGIS (csv files)
tide_dir = join(anuga_dir, 'tide_data')


# locations for boundary conditions
if source=='shark_bay':
    if boundary_event == '10000':
        boundary_file_name = 'shark_bay_3867_18052007'
        boundary_dir = join(anuga_dir, 'boundaries', 'shark_bay', '1_10000')
    elif boundary_event == 'july2006':
        boundary_file_name = 'july2006'    
        boundary_dir = join(anuga_dir, 'boundaries',
                            'shark_bay', 'july_2006_event')
    else:
        pass
        #raise Exception, 'Unknown boundary event specified'
        
if boundary_event != 'experimental':
    boundary_name = join(boundary_dir, boundary_file_name)
else:
    boundary_name = 'nil'

#output locations
output_dir = join(anuga_dir, 'outputs')+sep
output_build_time_dir = output_dir+build_time+sep
output_run_time_dir = output_dir +run_time+dir_comment+sep
output_run_time_dir_name = output_run_time_dir + scenario_name  #Used by post processing

#gauges
gauge_name = 'gaugesv3.csv' 
#gauge_name_test = 'gauge_checking_test.csv' 
gauges_dir = anuga_dir+sep+'gauges'+sep
gauges_dir_name = gauges_dir + gauge_name
#gauges_dir_name_test = gauges_dir + gauge_name_test

#tide_dir = anuga_dir+'tide_data'+sep

#buildings_filename = gauges_dir + 'pt_hedland_res.csv'
#buildings_filename_out = 'pt_hedland_res_modified.csv'
#buildings_filename_damage_out = 'pt_hedland_res_modified_damage.csv'
#tidal_filename = tide_dir + 'pt_hedland_tide.txt'
#community_filename = gauges_dir + 'CHINS_v2.csv'
#community_scenario = gauges_dir + 'community_pt_hedland.csv'

# clipping region to make DEM (pts file) from onshore data
#eastingmin = 594000
#eastingmax = 715000
#northingmin = 7720000
#northingmax = 7880000

# for ferret2sww
south = degminsec2decimal_degrees(-20,30,0)
north = degminsec2decimal_degrees(-17,10,0)
west = degminsec2decimal_degrees(117,00,0)
east = degminsec2decimal_degrees(120,00,0)

# region to export (used from export_results.py)
e_min_area = 713500   # Eastings min
e_max_area = 726350
n_min_area = 7100890
n_max_area = 7111250

#export_region = [[e_min_area,n_min_area],[e_min_area,n_max_area],
#                 [e_max_area,n_max_area],[e_max_area,n_min_area]]
                 
#FIXME (Ole): It is bad to read in polygons in project.py. If a file does not exist
#project.py cannot be imported

# Bounding polygon
bounding_polygon = read_polygon(join(polygons_dir, 'boundary_extent_small.csv'))
res_bounding_polygon = 500000*res_factor

# This is a restricted boundary used for experiments.
# Cannot use URS boundary on this!
small_bounding_polygon = read_polygon(join(polygons_dir, 'small_boundary.csv'))
res_bounding_polygon = 500000*res_factor


# Interior regions
interior_regions_and_resolutions = {
    'study_area': 300000,
    'bernier_dorre_islands': 250000,
    'map_area': 10000,
    'bay_area': 4000,
    'south_bank_buffer': 2000,
    'south_coast_polygon_1': 1000,
    'south_coast_polygon_2': 1000,
    'north_bay_coastline': 1000,
    'tsunami_approach_area': 50,
    'inundated_area': 10,
    'exclusion_1': 10000000,
    'exclusion_2': 10000000,
    'exclusion_3': 10000000}


#poly_exclusion_area1 = read_polygon(join(polygons_dir, 'exclusion_area1.csv'))
#res_exclusion_area1 = 2000000*res_factor
#
#interior_regions = [[poly_inundated_area, res_inundated_area],
#                    [poly_tsunami_approach_area, res_tsunami_approach_area],
#                    [poly_bay_area, res_bay_area],                    
#                    [poly_map_area, res_map_area]]
                    #[poly_model_area, res_model_area],
                    #[poly_exclusion_area1, res_exclusion_area1]]                                        

# FIXME (Ole): What if tags are wrong?

#Large bounding poly
#boundary_tags = {'back': [1, 2, 3, 4, 5],   
#                 'side': [0, 6],
#                 'ocean': [7, 8, 9, 10, 11]}


#Small bounding poly
#(counter clockwise starting in upper right corner)
boundary_tags = {'tide': [0, 1, 2],
                 'ocean': [3, 4, 5, 6]}

#Very small bounding poly (experimental boundary)
boundary_tags_small = {'tide': [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 12, 13],
                       'ocean': [11]}



# Read filenames etc (FIXME: move this to new file: read_structures)

interior_regions = [] # Consider using dictionary for mesh interface.
for area_name in interior_regions_and_resolutions.keys():
    file_name = join(polygons_dir, area_name + '.csv')
    
    print 'Reading polygon from', file_name 
    polygon = read_polygon(file_name)
    base_resolution = interior_regions_and_resolutions[area_name]
    
    interior_regions.append([polygon, base_resolution*res_factor])




trigs_min = number_mesh_triangles(interior_regions,
                                  bounding_polygon,
                                  res_bounding_polygon)



# For use with resetting IC (This is temporary)
poly_tsunami_approach_area = read_polygon(
    join(polygons_dir, 'tsunami_approach_area.csv'))



onshore_polygon = read_polygon(join(topographies_in_dir,
                                    'initial_condition.txt'))