source: anuga_work/production/patong/build_patong.py @ 6806

"""
Script for building the elevation data to run a tsunami inundation scenario
for patong, Thailand.

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

"""

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

# input elevation directory filenames
elevation_in_dir_name1 = project.elevation_in_dir_name1
elevation_in_dir_name2 = project.elevation_in_dir_name2
elevation_in_dir_name3 = project.elevation_in_dir_name3
elevation_in_dir_name4 = project.elevation_in_dir_name4

# output elevation directory filenames
elevation_dir_name1 = project.elevation_dir_name1
elevation_dir_name2 = project.elevation_dir_name2
elevation_dir_name3 = project.elevation_dir_name3
elevation_dir_name4 = project.elevation_dir_name4

# create elevation pts files
print'create Geospatial data1 objects from topographies'
G1 = Geospatial_data(file_name = elevation_in_dir_name1)
print'create Geospatial data2 objects from topographies'
G2 = Geospatial_data(file_name = elevation_in_dir_name2)
print'create Geospatial data3 objects from topographies'
G3 = Geospatial_data(file_name = elevation_in_dir_name3)
print'create Geospatial data4 objects from topographies'
G4 = Geospatial_data(file_name = elevation_in_dir_name4)

print 'add G1 + G2'
G5 = G1+G2

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

print 'clipping outside extents'
#G1_clip = G1.clip(project.extent_elev_dir4, verbose=True) 
G3_clip1 = G3.clip_outside(project.extent_elev_dir1, verbose=True) 
G3_clip2 = G3_clip1.clip_outside(project.extent_elev_dir2, verbose=True)
G4_clip = G4.clip_outside(project.extent_elev_dir3, verbose=True)

#print 'hello', dir(G3_clip1), dir(G3_clip1), dir(G4_clip)
print 'Add all geospatial objects'
G = G1 + G2 + G3_clip2 + G4_clip

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

print 'Oxport combined DEM file'
G_clipped.export_points_file(project.combined_dir_name + '.pts')
G_clipped.export_points_file(project.combined_dir_name + '.txt') #Use for comparision in ARC