"""Script for running a tsunami inundation scenario for Cairns, QLD Australia.

Source data such as elevation and boundary data is assumed to be available in
directories specified by project.py
The output sww file is stored in directory named after the scenario, i.e
slide or fixed_wave.

The scenario is defined by a triangular mesh created from project.polygon,
the elevation data and a tsunami wave generated by a submarine mass failure.

Ole Nielsen and Duncan Gray, GA - 2005 and Jane Sexton and
Nick Bartzis, GA - 2006
"""

def convert_arcgis_latlon_list_to_utm(points):
     #Used because arc gis produced csv files put lat lon in
     #reverse order to those accpeted by convert_latlon_to_utm()
     
     from anuga.coordinate_transforms.geo_reference import Geo_reference
     from anuga.coordinate_transforms.redfearn import redfearn
     old_geo = Geo_reference() 
     utm_points = []
     for point in points:
         zone, easting, northing = redfearn(float(point[1]),float(point[0]))
         new_geo = Geo_reference(zone)
         old_geo.reconcile_zones(new_geo) 
         utm_points.append([easting,northing])

     return utm_points, old_geo.get_zone()

#------------------------------------------------------------------------------
# Import necessary modules
#------------------------------------------------------------------------------

# Standard modules
import os
import time
import sys

# Related major packages
from anuga.shallow_water import Domain
from anuga.shallow_water import Reflective_boundary
from anuga.shallow_water import Dirichlet_boundary
from anuga.shallow_water import Time_boundary
from anuga.shallow_water import File_boundary
from anuga.pmesh.mesh_interface import create_mesh_from_regions
from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf, ferret2sww
from anuga.shallow_water.data_manager import dem2pts
from anuga.coordinate_transforms.redfearn import convert_from_latlon_to_utm
#from anuga.fit_interpolate.fit import fit_to_mesh_file (does not exist)

# Application specific imports
import project                 # Definition of file names and polygons


#------------------------------------------------------------------------------
# Define scenario as either slide or fixed_wave.
#------------------------------------------------------------------------------

#scenario = 'coseismic'
scenario = 'fixed_wave'

if os.access(scenario, os.F_OK) == 0:
    os.mkdir(scenario)
basename = scenario + 'source'


#------------------------------------------------------------------------------
# Preparation of topographic data
# Convert ASC 2 DEM 2 PTS using source data and store result in source data
#------------------------------------------------------------------------------

# Filenames
dem_name = 'boxingday' 

# Create DEM from asc data
#convert_dem_from_ascii2netcdf(dem_name, use_cache=True, verbose=True)
# Create pts file for onshore DEM
#dem2pts(dem_name, use_cache=True, verbose=True)


#------------------------------------------------------------------------------
# Create the triangular mesh based on overall clipping polygon with a tagged
# boundary and interior regions defined in project.py along with
# resolutions (maximal area of per triangle) for each polygon
#------------------------------------------------------------------------------

remainder_res = 10000000
islands_res = 100000
cairns_res = 100000
shallow_res = 500000
#interior_regions = [[project.poly_cairns, cairns_res]#,
#                    [project.poly_island0, islands_res],
#                    [project.poly_island1, islands_res],
#                    [project.poly_island2, islands_res],
#                    [project.poly_island3, islands_res],
#                    [project.poly_shallow, shallow_res]]

# filenames
meshname = project.meshname + '.tsh'
mesh_elevname = project.mesh_elevname  + '.tsh'

bounding_polygon, zone = convert_arcgis_latlon_list_to_utm(project.bounding_polygon_latlon)

from caching import cache
print 'start create mesh from regions'

if scenario == 'coseismic':
     _ = cache(create_mesh_from_regions,
               bounding_polygon,
               {'boundary_tags': {'one': [0], 'two': [1],
                                  'three': [2], 'four': [3],
                                  'five': [4], 'six': [5],
                                  'seven': [6], 'eight': [7],
                                  'nine': [8], 'ten': [9],
                                  'eleven': [10], 'twelve': [11],
                                  'thirteen': [12], 'fourteen': [13]},
                'maximum_triangle_area': 5000000,
                'filename': meshname},#,           
               #'interior_regions': interior_regions},
               verbose = True
               #, evaluate=True
               )

if scenario == 'fixed_wave':
     south = 7.40911081272
     north = 8.71484358635
     east = 99.1683687224
     west = 97.513856322
     points = [[south,west],[south,east],[north,east],[north,west]]

     bounding_polygon,zone=convert_from_latlon_to_utm(points)

     _ = cache(create_mesh_from_regions,
               bounding_polygon,
               {'boundary_tags': {'ocean_west': [0], 'bottom': [1],
               'onshore': [2], 'top': [3]},
                'maximum_triangle_area': 500000000,
                'filename': meshname},#,           
               #'interior_regions': interior_regions},
               verbose = True
               #, evaluate=True
               )

## cache(fit_to_mesh_file,(meshname,
##       project.combined_dir_name + '.pts',
##       mesh_elevname),
##       #{'verbose': True},
##       verbose = False
##      )

#------------------------------------------------------------------------------
# Setup computational domain
#------------------------------------------------------------------------------

domain = Domain(meshname, use_cache=False, verbose=True)

print 'Number of triangles = ', len(domain)
print 'The extent is ', domain.get_extent()
print domain.statistics()

domain.set_name(basename) 
domain.set_datadir(scenario)
domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
domain.set_minimum_storable_height(0.01)


#------------------------------------------------------------------------------
# Setup initial conditions
#------------------------------------------------------------------------------

def elevation(x,y):
     return -10.0
     
tide = 0.0
domain.set_quantity('stage', tide)
domain.set_quantity('friction', 0.0) 

if scenario == 'fixed_wave':
     # test with coarser bathymetry                        
     domain.set_quantity('elevation',
                         #filename = 'thaicoas_9.pts',
                         elevation,
                         use_cache=True,
                         verbose=True,
                         alpha=0.1)

if scenario == 'coseismic':
     domain.set_quantity('elevation', 
                         filename=project.combined_dir_name + '.pts',
                         use_cache=True,
                         verbose=True,
                         alpha=0.1)
#------------------------------------------------------------------------------
# Setup boundary conditions
#------------------------------------------------------------------------------

print 'start ferret2sww'
south = project.south 
north = project.north
west = project.west
east = project.east

"""
tide = 0.0
# Note only need to do when an SWW file for the MOST boundary doesn't exist
cache(ferret2sww,
      (project.boundary_most_in,
       project.boundary_most_out), 
      {'verbose': True,
       'minlat': south,
       'maxlat': north,
       'minlon': west,
       'maxlon': east,
#       'origin': project.mesh_origin,
       'origin': domain.geo_reference.get_origin(),
       'mean_stage': tide,
       'zscale': 10,                 #Enhance tsunami (Does this affect result)
       'fail_on_NaN': False,
       'inverted_bathymetry': True},
       #evaluate = True,
       verbose = True)#,
      #dependencies =  'most_results' + '.sww')

"""
print 'Available boundary tags', domain.get_boundary_tags()

#Bf = File_boundary('most_results'+'.sww', 
#                    domain, verbose = True)
Br = Reflective_boundary(domain)
Bd = Dirichlet_boundary([tide,0,0])
Bw = Time_boundary(domain = domain,
                       f=lambda t: [(60<t<3660)*50, 0, 0])
    
# 60 min square wave starting at 1 min, 50m high
if scenario == 'fixed_wave':
    domain.set_boundary({'ocean_west': Br,
                         'bottom': Br,
                         'onshore': Br,
                         'top': Br})

# boundary conditions for slide scenario
if scenario == 'coseismic':
    domain.set_boundary( {'one': Bf,'two': Bf, 'three':  Bf, 'four': Bf,
                      'five': Bf, 'six': Bd, 'seven': Bd, 'eight': Bf,
                      'nine': Bf, 'ten':  Bd, 'eleven': Bd, 'twelve': Bf,
                      'thireteen': Bf, 'fourteen': Bf})
    

#------------------------------------------------------------------------------
# Evolve system through time
#------------------------------------------------------------------------------
#domain.visualise = True
#domain.visualise_color_stage = True

import time
t0 = time.time()

from Numeric import allclose
from anuga.abstract_2d_finite_volumes.quantity import Quantity

if scenario == 'cosesimic':
    
     for t in domain.evolve(yieldstep = 10, finaltime = 60):
          print domain.quantities['elevation'].vertex_values
          domain.write_time() 

     for t in domain.evolve(yieldstep = 600, finaltime = 21600, 
                            skip_initial_step = True):
        domain.write_time()
        domain.write_boundary_statistics(tags = 'ocean_west')

if scenario == 'fixed_wave':

    # save every two mins leading up to wave approaching land
    for t in domain.evolve(yieldstep = 10, finaltime = 50):
         print domain.quantities['elevation'].vertex_values
         print domain.quantities['stage'].vertex_values
         
         domain.write_time()

    #save every 30 secs as wave starts inundating ashore
    #for t in domain.evolve(yieldstep = 10, finaltime = 10000, 
    #                       skip_initial_step = True):
    #    domain.write_time()
            
print 'That took %.2f seconds' %(time.time()-t0)