source: anuga_work/production/perth/build_perth_250m.py @ 5782

"""
Script for building the elevation data to run tsunami inundation scenario
for Perth, WA, Australia.

Input: elevation data from project_250m.py
Output: pts file stored in project_250m.topographies_dir 
The run_perth.py is reliant on the output of this script.

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.data_manager import convert_dem_from_ascii2netcdf, dem2pts
from anuga.geospatial_data.geospatial_data import *
from anuga.shallow_water.data_manager import start_screen_catcher, copy_code_files

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

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

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

start_screen_catcher(project_250m.output_build_time_dir)

print 'USER:    ', project_250m.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_250m.poly_all",project_250m.poly_all
print"project_250m.combined_dir_name",project_250m.combined_dir_name

# input elevation directory filenames
onshore_in_dir_name = project_250m.onshore_in_dir_name

# output elevation directory filenames
onshore_dir_name = project_250m.onshore_dir_name

# creates DEM from asc data
print "creates DEMs from asc data"
convert_dem_from_ascii2netcdf(onshore_in_dir_name, basename_out=onshore_dir_name, use_cache=True, verbose=True)

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

#-------------------------------------------------------------------------------
# Combine datasets into project_250m.combined_dir_name
#-------------------------------------------------------------------------------

print'create Geospatial data1 objects from topographies'
G = Geospatial_data(file_name = onshore_dir_name + '.pts')

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

print'export combined DEM file'
if access(project_250m.topographies_dir,F_OK) == 0:
    mkdir (project_250m.topographies_dir)
G_clipped.export_points_file(project_250m.combined_dir_name + '.pts')
#G_clipped.export_points_file(project_250m.combined_dir_name + '.txt') #Use for comparision in ARC

import sys
sys.exit()