source: anuga_work/production/broome/build_broome.py @ 5001

"""Script for running tsunami inundation scenario for Dampier, WA, Australia.

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

The scenario is defined by a triangular mesh created from project_urs.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.shallow_water.data_manager import start_screen_catcher, copy_code_files
from anuga_parallel.parallel_abstraction import get_processor_name

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

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

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

start_screen_catcher(project_urs.output_build_time_dir)

print 'USER:    ', project_urs.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_urs.combined_dir_name",project_urs.combined_dir_name

# topography directory filenames
onshore_in_dir_name = project_urs.onshore_in_dir_name
onshore_in_dir_name1 = project_urs.onshore_in_dir_name1
coast_in_dir_name = project_urs.coast_in_dir_name
#island_in_dir_name = project_urs.island_in_dir_name
offshore_in_dir_name = project_urs.offshore_in_dir_name
offshore_in_dir_name1 = project_urs.offshore_in_dir_name1
offshore_in_dir_name2 = project_urs.offshore_in_dir_name2
offshore_in_dir_name3 = project_urs.offshore_in_dir_name3

onshore_dir_name = project_urs.onshore_dir_name
onshore_dir_name1 = project_urs.onshore_dir_name1
coast_dir_name = project_urs.coast_dir_name
#island_dir_name = project_urs.island_dir_name
offshore_dir_name = project_urs.offshore_dir_name
offshore_dir_name1 = project_urs.offshore_dir_name1
offshore_dir_name2 = project_urs.offshore_dir_name2
offshore_dir_name3 = project_urs.offshore_dir_name3

# 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(onshore_in_dir_name1, basename_out=onshore_dir_name1, use_cache=True, verbose=True)
convert_dem_from_ascii2netcdf(offshore_in_dir_name1, basename_out=offshore_dir_name1, use_cache=True, verbose=True)
convert_dem_from_ascii2netcdf(offshore_in_dir_name2, basename_out=offshore_dir_name2, use_cache=True, verbose=True)
convert_dem_from_ascii2netcdf(offshore_in_dir_name3, basename_out=offshore_dir_name3, 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(onshore_dir_name1, use_cache=True, verbose=True)
#creates pts file for island DEM
dem2pts(offshore_dir_name1, use_cache=True, verbose=True)
dem2pts(offshore_dir_name2, use_cache=True, verbose=True)
dem2pts(offshore_dir_name3, use_cache=True, verbose=True)

print'create Geospatial data1 objects from topographies',onshore_dir_name + '.pts'
G1 = Geospatial_data(file_name = onshore_dir_name + '.pts')
print'create Geospatial data1a objects from topographies',onshore_dir_name1 + '.pts'
G1a = Geospatial_data(file_name = onshore_dir_name1 + '.pts')
print'create Geospatial data2 objects from coast', coast_in_dir_name + '.txt'
G2 = Geospatial_data(file_name = coast_in_dir_name + '.txt')
#print'create Geospatial data3 objects from island'
#G3 = Geospatial_data(file_name = island_dir_name + '.pts')
print'create Geospatial data3 objects from offshore',offshore_in_dir_name + '.txt'
G_off = Geospatial_data(file_name = offshore_in_dir_name + '.txt')
print'create Geospatial data4 objects from offshore1',offshore_dir_name1 + '.pts'
G_off1 = Geospatial_data(file_name = offshore_dir_name1 + '.pts')
print'create Geospatial data4 objects from offshore2',offshore_dir_name2 + '.pts'
G_off2 = Geospatial_data(file_name = offshore_dir_name2 + '.pts')
print'create Geospatial data objects from offshore3',offshore_dir_name3 + '.pts'
G_off3 = Geospatial_data(file_name = offshore_dir_name3 + '.pts')

print'add all geospatial objects'
G = G1 + G1a + G2 + G_off + G_off1 + G_off2 + G_off3

print'clip combined geospatial object by bounding polygon'
#G_clipped = G.clip(project_urs.poly_all)
#FIXME: add a clip function to pts
#print'shape of clipped data', G_clipped.get_data_points().shape

print'export combined DEM file'
if access(project_urs.topographies_dir,F_OK) == 0:
    mkdir (project_urs.topographies_dir)
print'export',project_urs.combined_dir_name+ '.txt'
G.export_points_file(project_urs.combined_dir_name+ '.txt')

'''
print'split'
G_small, G_other = G.split(.10,verbose=True)


print 'export',project_urs.combined_dir_name + '.txt'
G.export_points_file(project_urs.combined_dir_name + '.txt')
print 'export', project_urs.combined_small_dir_name + '.txt'
G_small.export_points_file(project_urs.combined_small_dir_name + '.txt')
#G_clipped.export_points_file(project_urs.combined_dir_name + '.xya')



print'export',project_urs.combined_dir_name+ '.txt'
G_all=Geospatial_data(file_name = project_urs.combined_dir_name+ '.txt')
print'split'
G_all_1, G_all_2 = G_all.split(.10)
print'export 1'
G_all_1.export_points_file(project_urs.combined_dir_name+'_small' + '.xya')
print'export 2'
G_all_2.export_points_file(project_urs.combined_dir_name+'_other1' + '.xya')


#-------------------------------------------------------------------------
# Convert URS to SWW file for boundary conditions 
#-------------------------------------------------------------------------
print 'starting to create boundary conditions'
from anuga.shallow_water.data_manager import urs2sww, urs_ungridded2sww

print 'boundaries_in_dir_name',project_urs.boundaries_in_dir_name

urs_ungridded2sww(project_urs.boundaries_in_dir_name, project_urs.boundaries_in_dir_name,
                  verbose=True, mint=4000, maxt=35000, zscale=1)


'''