source: anuga_work/production/carnarvon/build_carnarvon.py @ 6514

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

Input: elevation data from project.py
Output: pts file stored in project.topographies_dir 
The run_geraldton.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
onshore_in_dir_name = project.onshore_in_dir_name
offshore_in_dir_name = project.offshore_in_dir_name
offshore_in_dir_name1 = project.offshore_in_dir_name1
offshore_in_dir_aoi = project.offshore_in_dir_aoi
offshore_in_dir_aos = project.offshore_in_dir_aos

# output elevation directory filenames
onshore_dir_name = project.onshore_dir_name
offshore_dir_name = project.offshore_dir_name
offshore_dir_name1 = project.offshore_dir_name1
offshore_dir_aoi = project.offshore_dir_aoi
offshore_dir_aos = project.offshore_dir_aos


# 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)
convert_dem_from_ascii2netcdf(offshore_in_dir_aoi, basename_out=offshore_dir_aoi, use_cache=True, verbose=True)
convert_dem_from_ascii2netcdf(offshore_in_dir_aos, basename_out=offshore_dir_aos, 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)
dem2pts(offshore_dir_aoi ,use_cache=True,verbose=True)
dem2pts(offshore_dir_aos ,use_cache=True,verbose=True)



# create onshore pts files
print'create Geospatial onshore objects from topographies'
G1 = Geospatial_data(file_name = onshore_dir_name + '.pts')

# create coastal and offshore txt files
print'create Geospatial offshore objects from topographies'
G_off1 = Geospatial_data(file_name = offshore_in_dir_name)
print'create Geospatial makeup objects from topographies'
G_off2 = Geospatial_data(file_name = offshore_in_dir_name1)
print'create Geospatial aoi objects from topographies'
G_offaoi = Geospatial_data(file_name = offshore_dir_aoi+ '.pts')
print'create Geospatial aos objects from topographies'
G_offaos = Geospatial_data(file_name = offshore_dir_aos+ '.pts')


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

print 'add all bathymetry files'
G_off = G_off1 + G_off2

##print 'clipping data to coast'
##G_ocean = G_off.clip(project.poly_ocean, verbose=True) #mainland
##G_ocean_b = G_ocean.clip_outside(project.poly_island1, verbose=True) #bernier
##G_ocean_bd = G_ocean_b.clip_outside(project.poly_island2, verbose=True) #dorne
##
##G1_clip_m = G1.clip(project.poly_mainland, verbose=True)
##G1_clip_b = G1.clip(project.poly_island1, verbose=True)
##G1_clip_d = G1.clip(project.poly_island2, verbose=True)
##G1_clip = G1_clip_m + G1_clip_b + G1_clip_d

print'add all geospatial objects'
G_aos = G1 + G_off
G_aosclip = G_aos.clip_outside(project.poly_aos_extract, verbose=True)
G_aoi = G_aosclip.clip_outside(project.poly_aoi1, verbose=True)
G = G_aoi +G_offaoi+G_offaos

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