source: anuga_work/production/onslow_2006/build_onslow_grad.py @ 6855

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

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

The scenario is defined by a triangular mesh created from project_grad_urs.polygon,
the elevation data and a simulated tsunamigenic earthquake.

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_grad   # Definition of file names and polygons


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

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

start_screen_catcher(project_grad.output_build_time_dir)

print 'Processor Name:', get_processor_name()
print 'User: ', project_grad.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_grad.combined_dir_name', project_grad.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_grad.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_grad.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_grad.bounding_polygon)


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

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

print 'boundary_dir', project_grad.boundaries_in_dir

maxt = 40000
urs_ungridded2sww(project_grad.boundaries_in_dir+sep+project_grad.boundaries_name,
                  project_grad.boundaries_in_dir+sep+project_grad.boundaries_name+'_'+ str(maxt),
                  verbose=True,
                  mint=0, maxt=maxt, zscale=1)

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