source: anuga_work/production/bunbury_storm_surge_2009/build_elevation.py @ 7629

"""Build the elevation data to run a storm surge inundation scenario 
for Mandurah, Western Australia, 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)

        G_grid = Geospatial_data(file_name=absolute_filename+'.pts',
                                                    verbose=True)
        print 'Clip geospatial object'
        geospatial_data[filename] = G_grid.clip(project.bounding_polygon)

# Create Geospatial data from TXT files
if not project.point_filenames == []:
    for filename in project.point_filenames:
        absolute_filename = join(project.topographies_folder, filename)
        G_points = Geospatial_data(file_name=absolute_filename,
                                                    verbose=True)
        print 'Clip geospatial object'
        geospatial_data[filename] = G_points.clip(project.bounding_polygon)

print 'Geospatial data objects: ', geospatial_data.keys()

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

print 'Add geospatial objects' # except', project.offshore_name5

print 'project.elevation_clip_box', project.elevation_clip_box

keylist = project.ascii_grid_filenames + project.point_filenames

G1 = geospatial_data['swwa_10m'].clip_outside(project.elevation_clip_box)
G2 = geospatial_data['bunbury_5m'].clip(project.elevation_clip_box)
#~ G2 = geospatial_data['bunbury_nth_a'].clip(project.elevation_clip_box)
#~ G3 = geospatial_data['bunbury_nth_b'].clip(project.elevation_clip_box)
#~ G4 = geospatial_data['bunbury_sth_a'].clip(project.elevation_clip_box)
#~ G5 = geospatial_data['bunbury_sth_b'].clip(project.elevation_clip_box)
G3 = geospatial_data['DPI.txt']
G4 = geospatial_data['Busselton_Chart_Clip_ss.txt']
G5 = geospatial_data['Busselton_NavyFinal_Clip_ss.txt']
G6 = geospatial_data['Leschenault_TIN.txt']
G7 = geospatial_data['DPI5U1A02_01a_edited.txt']
G8 = geospatial_data['DPI5U1A02_01b_edited.txt']
G9 = geospatial_data['DPI5U1A02_01c_edited.txt']
G10 = geospatial_data['DPI5U1A02_01d_edited.txt']
G11 = geospatial_data['DPI5U1A02_01e_edited.txt']

#~ G_list = [G1, G2, G3, G4, G5, G6, G7, G8, G9, G10]

#~ print 'G1', G1.attributes
#~ print 'G2', G2.attributes
#~ print 'G3', G3.attributes
#~ print 'G4', G4.attributes
#~ print 'G5', G5.attributes
#~ print 'G6', G6.attributes
#~ print 'G7', G7.attributes
#~ print 'G8', G8.attributes
#~ print 'G9', G9.attributes
#~ print 'G10', G10.attributes

G = G1 + G2 + G3 + G4 + G5 + G6 + G7 + G8 + G9 + G10 + G11

# G = None
# for key in geospatial_data:
    # if key == project.ascii_grid_filenames[0]:
                # G == geospatial_data[key].clip_outside(project.elevation_clip_box)
        # elif key == project.ascii_grid_filenames[1]:
                # G = geospatial_data[key].clip(project.elevation_clip_box)
    # continue
    
    # G += geospatial_data[key]

print 'Export combined DEM file'
G.export_points_file(project.combined_elevation + '.pts')

print 'Do txt version too'
# Use for comparision in ARC
#~ for g in G_list:
        #~ for i in keylist:
        #~ g.export_points_file(join(project.topographies_folder, (str(keylist[i]) +'_export.txt')))
G.export_points_file(project.combined_elevation + '.txt')
G1.export_points_file(join(project.topographies_folder, keylist[0] +'_export.txt'))
G2.export_points_file(join(project.topographies_folder, keylist[1] +'_export.txt'))
G3.export_points_file(join(project.topographies_folder, keylist[2] +'_export.txt'))
G4.export_points_file(join(project.topographies_folder, keylist[3] +'_export.txt'))
G5.export_points_file(join(project.topographies_folder, keylist[4] +'_export.txt'))
G6.export_points_file(join(project.topographies_folder, keylist[5] +'_export.txt'))
G7.export_points_file(join(project.topographies_folder, keylist[6] +'_export.txt'))
G8.export_points_file(join(project.topographies_folder, keylist[7] +'_export.txt'))
G9.export_points_file(join(project.topographies_folder, keylist[8] +'_export.txt'))
G10.export_points_file(join(project.topographies_folder, keylist[9] +'_export.txt'))
G11.export_points_file(join(project.topographies_folder, keylist[10] +'_export.txt'))