source: anuga_work/production/galle_2007/build_gal5.py @ 6806

"""

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

The scenario is defined by a triangular mesh created from project.polygon,
the elevation data and a simulated submarine landslide.

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.abstract_2d_finite_volumes.util 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
#-------------------------------------------------------------------------------
# Create DEM from asc data
'''
convert_dem_from_ascii2netcdf(galle_dem_name, use_cache=True, verbose=True)
convert_dem_from_ascii2netcdf(hikka_dem_name, use_cache=True, verbose=True)
convert_dem_from_ascii2netcdf(topo_dem_name, use_cache=True, verbose=True)
convert_dem_from_ascii2netcdf(bathycoa_dem_name, use_cache=True, verbose=True)
#convert_dem_from_ascii2netcdf(bathyfine_dem_name, use_cache=True, verbose=True)

# Create pts file for onshore DEM
dem2pts(galle_dem_name, use_cache=True, verbose=True)
dem2pts(hikka_dem_name, use_cache=True, verbose=True)
dem2pts(topo_dem_name, use_cache=True, verbose=True)
dem2pts(bathycoa_dem_name, use_cache=True, verbose=True)
#dem2pts(bathyfine_dem_name, use_cache=True, verbose=True)

# DEM
G1 = Geospatial_data(file_name = project.galle_dem_name + '.pts')
print 'loading fine bathy data'
G2 = Geospatial_data(file_name = project.hikka_dem_name + '.pts')
print 'loading fine bathy data'
G3 = Geospatial_data(file_name = project.topo_dem_name + '.pts')
G6 = Geospatial_data(file_name = project.coast_name + '.txt')
print 'loading fine bathy data'
G4 = Geospatial_data(file_name = project.bathyfine_dem_name + '.pts')
G5 = Geospatial_data(file_name = project.bathycoa_dem_name + '.pts')
G7 = Geospatial_data(file_name = project.canal_dem_name + '.txt')
G8 = Geospatial_data(file_name = project.fort_dem_name+'.txt')

# coastline
#G4.export_points_file(project.bathyfine_dem_name + '.txt')

print 'clip G5'
G5_clip = G5.clip_outside(Geospatial_data(project.poly_clip))
print 'adding data files together'
#G = G1 + G2 + G3 + G4 + G6 + G5_clip
G = G1 + G2 + G3 + G4 + G6 + G5 + G7 + G8
#G = G1 + G2 + G3 + G4 + G6 + G5 + G7
#G = G1 + G2 + G3 + G5 + G6
print 'exporting elevation file'
G.export_points_file(project.combined_canal_fort_dem_name + '.txt')


print'load data'
G = Geospatial_data(file_name = project.coast_name + '.txt')
#G = Geospatial_data(file_name = project.combined_canal_dem_name + '.txt')
print'split combined data set'
G_small, G_other = G.split(0.1, True)
G_small.export_points_file(project.combined_small_dem_name+ '1.txt')
'''
# if you would like to look at your points file in GIS environment
# use this line and then change extension to txt for import into ArcGIS
#G.export_points_file(combined_dem_name + '.xya')


print 'start ferret2sww'
from anuga.shallow_water.data_manager import ferret2sww

#note only need to do when an SWW file for the MOST boundary doesn't exist

ferret2sww(project.boundaries_in_dir_name,
       project.boundaries_dir_name,
       verbose = True,
#       minlat = project.south_source, 
#       maxlat = project.north_source,
#       minlon = project.west_source, 
#       maxlon = project.east_source,
       mint = None, maxt = None,
       zscale = 1,
       fail_on_NaN = False,
       )