source: anuga_work/production/australia_ph2/ceduna/build_elevation.py @ 6476

"""Build the elevation data to run a tsunami inundation scenario 
for Busselton, WA, Australia.

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

"""

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

# Standard modules
import os
from os.path import join

# Related major packages
from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf
from anuga.shallow_water.data_manager import dem2pts
from anuga.shallow_water.data_manager import start_screen_catcher
from anuga.shallow_water.data_manager import copy_code_files
from anuga.geospatial_data.geospatial_data import Geospatial_data

# Application specific imports
from setup_model 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,__file__, 
                os.path.dirname(project.__file__)+os.sep+\
                project.__name__+'.py' )
start_screen_catcher(project.output_build)


#------------------------------------------------------------------------------
# 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', project.bounding_polygon
print 'project.combined_elevation_basename', project.combined_elevation_basename

# Create Geospatial data from ASCII files
geospatial_data = {}
if not project.ascii_grid_filenames == []:
    for filename in project.ascii_grid_filenames:
        absolute_filename = join(project.topographies_folder, filename)
        convert_dem_from_ascii2netcdf(absolute_filename,
                                      basename_out=absolute_filename,
                                      use_cache=True,
                                      verbose=True)
        dem2pts(absolute_filename, use_cache=True, verbose=True)

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

# Create Geospatial data from TXT files
if not project.point_filenames == []:
    for filename in project.point_filenames:
        absolute_filename = join(project.topographies_folder, filename)
        geospatial_data[filename] = Geospatial_data(file_name=absolute_filename,
                                                    verbose=True)

#-------------------------------------------------------------------------------
# Combine, clip and export dataset 
#-------------------------------------------------------------------------------

print 'Add geospatial objects' # except', project.offshore_name5
G = None
for key in geospatial_data:
    #if key == project.point_filenames[4] or key == project.ascii_filenames[1]:
    #    # Skip these files for now
    #    continue
    
    G += geospatial_data[key]
        
print 'Clip combined geospatial data'
##G_clip = G.clip_outside(project.poly_aoi1)
##G_all = G_clip + geospatial_data[project.point_filenames[4]] 
#G_clipped = G_all.clip(project.poly_all)
G_clip = G.clip(project.bounding_polygon)


print 'Export combined DEM file'
G_clip.export_points_file(project.combined_elevation + '.pts')
print 'Do txt version too'
# Use for comparision in ARC
G_clip.export_points_file(project.combined_elevation + '.txt')