source: anuga_work/development/boxingday_2007/build_boxing.py @ 4297

"""Script for running tsunami inundation scenario for boxing day Phuket,
Thailand. 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 project.output_time_dir 

The scenario is defined by a triangular mesh created from project.polygon,
the elevation data and a simulated submarine landslide.

John Jakeman - ANU 2007
"""

#use _read_asc() in data_manager


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

# Standard modules
from os import sep
from os.path import dirname, basename
from os import mkdir, access, F_OK
from shutil import copy
import time
import sys
from time import localtime, strftime

# Related major packages
from anuga.shallow_water import Domain
from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf, dem2pts, ferret2sww
from anuga.geospatial_data.geospatial_data import *
from anuga.abstract_2d_finite_volumes.util import start_screen_catcher, copy_code_files
from anuga.utilities.polygon import polygon_area

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

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()
"""
#------------------------------------------------------------------------------
# Preparation of topographic data
#
# Convert ASC 2 DEM 2 PTS using source data and store result in source data
# Do for coarse and fine data
# Fine pts file to be clipped to area of interest
#------------------------------------------------------------------------------
print"project.bounding_polygon_latlon",project.bounding_polygon_latlon

# topography and bathymetry directory filenames
bathymetry_file = project.bathymetry

#gets time for new dir
time = strftime('%Y%m%d_%H%M%S',localtime())

# creates DEM from asc data
print "creates DEMs from asc data"
convert_dem_from_ascii2netcdf(bathymetry_file, basename_out=project.bathymetry_points, use_cache=True, verbose=True)

#creates pts file for onshore DEM
print "creates pts file for onshore DEM"
dem2pts(project.bathymetry_points,
        use_cache=True, 
        verbose=True)

print'create Geospatial data objects from topographies'
G = Geospatial_data(file_name = project.bathymetry_points + '.pts')

bounding_polygon, zone = convert_arcgis_latlon_list_to_utm(project.bounding_polygon_latlon)
print bounding_polygon
print 'zone of bounding polygon', zone
print 'Area of bounding polygon', polygon_area(bounding_polygon)/1000000.0

print'clip combined geospatial object by bounding polygon'
G_clipped = G.clip(bounding_polygon)

print'export combined DEM file'
G_clipped.export_points_file(project.combined_dir_name + '.pts')
"""