Index: /anuga_work/production/patong/build_patong.py
===================================================================
--- /anuga_work/production/patong/build_patong.py	(revision 6269)
+++ /anuga_work/production/patong/build_patong.py	(revision 6269)
@@ -0,0 +1,101 @@
+"""
+Script for building the elevation data to run a tsunami inundation scenario
+for patong, Thailand.
+
+Input: elevation data from project.py
+Output: pts file stored in project.topographies_dir 
+The run_patong.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
+elevation_in_dir_name1 = project.elevation_in_dir_name1
+elevation_in_dir_name2 = project.elevation_in_dir_name2
+elevation_in_dir_name3 = project.elevation_in_dir_name3
+elevation_in_dir_name4 = project.elevation_in_dir_name4
+
+# output elevation directory filenames
+elevation_dir_name1 = project.elevation_dir_name1
+elevation_dir_name2 = project.elevation_dir_name2
+elevation_dir_name3 = project.elevation_dir_name3
+elevation_dir_name4 = project.elevation_dir_name4
+
+# create elevation pts files
+print'create Geospatial data1 objects from topographies'
+G1 = Geospatial_data(file_name = elevation_in_dir_name1)
+print'create Geospatial data2 objects from topographies'
+G2 = Geospatial_data(file_name = elevation_in_dir_name2)
+print'create Geospatial data3 objects from topographies'
+G3 = Geospatial_data(file_name = elevation_in_dir_name3)
+print'create Geospatial data4 objects from topographies'
+G4 = Geospatial_data(file_name = elevation_in_dir_name4)
+
+print 'add G1 + G2'
+G5 = G1+G2
+
+#-------------------------------------------------------------------------------
+# Combine, clip and export dataset 
+#-------------------------------------------------------------------------------
+
+print 'clipping outside extents'
+G1_clip = G1.clip(project.extent_elev_dir4, verbose=True) 
+G3_clip1 = G3.clip_outside(project.extent_elev_dir4, verbose=True) 
+G3_clip2 = G3_clip1.clip_outside(project.extent_elev_dir2, verbose=True)
+G4_clip = G4.clip_outside(project.extent_elev_dir3, verbose=True)
+
+#print 'hello', dir(G3_clip1), dir(G3_clip1), dir(G4_clip)
+print'add all geospatial objects'
+G = G1_clip + G2 + G3_clip1 + G3_clip2 + G4_clip
+
+print'clip combined geospatial object by bounding polygon'
+G_clipped = G.clip(project.poly_all)
+
+print'export combined DEM file'
+if access(project.topographies_dir,F_OK) == 0:
+    mkdir (project.topographies_dir)
+G_clipped.export_points_file(project.combined_dir_name + '.pts')
+G_clipped.export_points_file(project.combined_dir_name + '.txt') #Use for comparision in ARC
+
