"""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

# G1 = geospatial_data['Man_25m'].clip(project.elevation_clip_box)
# G2 = geospatial_data['ph10m_ss'].clip_outside(project.elevation_clip_box)

G1 = geospatial_data['m_peel_aoi'].clip(project.elevation_clip_box)
G2 = geospatial_data['ph10m_ss'].clip_outside(project.elevation_clip_box)
G3 = geospatial_data['DpiU1A03_ss.csv']
G4 = geospatial_data['MA-46893-SNDS_AHD.csv']
G5 = geospatial_data['MS0205HY_AHD.csv']
G6 = geospatial_data['MS0404_AHD.csv']
G7 = geospatial_data['YU0403HY_AHD.csv']
G8 = geospatial_data['original_data_ss.csv']

G = G1 + G2 + G3 + G4 + G5 + G6 + G7 + G8
	
print 'Export combined DEM file'
G.export_points_file(project.combined_elevation + '.pts')

print 'Do txt version too'
# Use for comparision in ARC
G.export_points_file(project.combined_elevation + '.txt')
G1.export_points_file(join(project.topographies_folder, (str(project.ascii_grid_filenames[0]) +'_export.txt')))
G2.export_points_file(join(project.topographies_folder, (str(project.ascii_grid_filenames[1]) + '_export.txt')))