source: anuga_work/production/shark_bay_2007/build_shark_bay.py @ 5563

"""Script for running tsunami inundation scenario for Dampier, WA, Australia.

Source data such as elevation and boundary data is assumed to be available in
directories specified by project_urs.py
The output sww file is stored in project.output_time_dir 

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

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

#------------------------------------------------------------------------------
# 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

# Related major packages
from anuga.shallow_water import Domain
from anuga.shallow_water import Dirichlet_boundary
from anuga.shallow_water import File_boundary
from anuga.shallow_water import Reflective_boundary
from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf
from anuga.shallow_water.data_manager import dem2pts
from anuga.pmesh.mesh_interface import create_mesh_from_regions
from anuga.geospatial_data.geospatial_data import Geospatial_data
from anuga.shallow_water.data_manager import start_screen_catcher
from anuga.shallow_water.data_manager import copy_code_files
from anuga.shallow_water.data_manager import urs_ungridded2sww
from anuga_parallel.parallel_abstraction import get_processor_name

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


#------------------------------------------------------------------------------
# Copy scripts to time stamped output directory and capture screen
# output to file
#------------------------------------------------------------------------------

copy_code_files(project.output_build_time_dir,__file__, 
               dirname(project.__file__)+sep+ project.__name__+'.py' )#

start_screen_catcher(project.output_build_time_dir)

print 'Processor Name:', get_processor_name()
print 'User: ', project.user


#-------------------------------------------------------------------------------
# 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.combined_dir_name', project.combined_dir_name


geospatial_data = None
# create DEMs from asc data

# FIXME (Ole): Clip data by interior regions if possible

print 'creating geospatial data objects from asc data (via dem and pts formats)'
for filename in project.ascii_grid_filenames:
    convert_dem_from_ascii2netcdf(filename,
                                  basename_out=filename,
                                  use_cache=True, verbose=True)
    dem2pts(filename, use_cache=True, verbose=True)

    geospatial_data += Geospatial_data(file_name=filename + '.pts',
                                       verbose=True)


print 'creating geospatial data objects from txt data'
for filename in project.point_filenames:
    geospatial_data += Geospatial_data(file_name=filename + '.txt',
                                       verbose=True)


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


print 'export combined geospatial data'
if access(project.topographies_dir, F_OK) == 0:
    mkdir (project.topographies_dir)
G.export_points_file(project.combined_dir_name + '.pts')


#print 'split combined data set'
#G_small, _ = G.split(0.1, True)
#
#print 'export sub sampled DEM file'
#G_small.export_points_file(project.combined_dir_name + '_small' + '.txt')


 
#-------------------------------------------------------------------------
# Convert URS to SWW file for boundary conditions 
#-------------------------------------------------------------------------
print 'converting boundary conditions to sww format'

print 'boundary_dir', project.boundary_dir
print 'project.boundary_name', project.boundary_name

urs_ungridded2sww(project.boundary_name, verbose=True,
                  mint=0, maxt=40000, zscale=1)


print 'Finished building the %s scenario' %project.scenario_name