Index: /anuga_work/production/busselton/busselton_rerun/Arc_txt2feature_class.py
===================================================================
--- /anuga_work/production/busselton/busselton_rerun/Arc_txt2feature_class.py	(revision 7118)
+++ /anuga_work/production/busselton/busselton_rerun/Arc_txt2feature_class.py	(revision 7119)
@@ -26,5 +26,5 @@
 time_dir6 = '20090330_120458_run_final_0_27283_Bt_kvanputt'
 
-time_dirs = [time_dir6] #1, time_dir2] #, time_dir5] #4, time_dir5, time_dir6]
+time_dirs = [time_dir1, time_dir4] #, time_dir5] #4, time_dir5, time_dir6]
 
 for time_dir in time_dirs:
Index: /anuga_work/production/busselton/busselton_rerun/Busselton4distribution.py
===================================================================
--- /anuga_work/production/busselton/busselton_rerun/Busselton4distribution.py	(revision 7119)
+++ /anuga_work/production/busselton/busselton_rerun/Busselton4distribution.py	(revision 7119)
@@ -0,0 +1,140 @@
+# ---------------------------------------------------------------------------
+# Carnarvon4distribution.py
+# Created on: Thu Apr 30 2009 03:51:03 PM
+#   (generated by ArcGIS/ModelBuilder)
+# ---------------------------------------------------------------------------
+
+# Import system modules
+import sys, string, os, arcgisscripting
+
+# Create the Geoprocessor object
+gp = arcgisscripting.create()
+
+# Check out any necessary licenses
+gp.CheckOutExtension("spatial")
+
+# Load required toolboxes...
+gp.AddToolbox("C:/Program Files/ArcGIS/ArcToolbox/Toolboxes/Spatial Analyst Tools.tbx")
+gp.AddToolbox("C:/Program Files/ArcGIS/ArcToolbox/Toolboxes/Data Management Tools.tbx")
+gp.overwriteoutput = 1
+
+community = 'Busselton'
+output_loc = "\\\\nas2\\gemd\\georisk_models\\inundation\\data\\western_australia\\FESA_outputs_Dec_08\\ARCgrids\\April_2009\\"
+input_scenario = "\\\\nas2\\gemd\\georisk_models\\inundation\\data\\western_australia\\busselton_tsunami_scenario\\"
+input_output = input_scenario + "\\anuga\\outputs\\"
+
+locations = ['Bunbury', 'Busselton']
+
+in_Tsunami_1_HAT = input_output + "20090327_152717_run_final_0.6_27283_Bt_kvanputt\\raster.gdb\\"
+in_Tsunami_1_MSL = input_output + "20090330_120458_run_final_0_27283_Bt_kvanputt\\raster.gdb\\"
+in_Tsunami_2_HAT = input_output + "20090324_134328_run_final_0.6_27255_Bt_kvanputt\\raster.gdb\\"
+in_Tsunami_2_MSL = input_output + "20090327_160646_run_final_0_27255_Bt_kvanputt\\raster.gdb\\"
+in_Tsunami_3_HAT = input_output + "20090324_140747_run_final_0.6_68693_Bt_kvanputt\\raster.gdb\\"
+in_Tsunami_3_MSL = input_output + "20090327_161520_run_final_0_68693_Bt_kvanputt\\raster.gdb\\"
+
+Tsunami_1_HAT = output_loc + "Tsunami1-Exmouth\\" + community + "_HAT.gdb\\"
+Tsunami_1_MSL = output_loc + "Tsunami1-Exmouth\\" + community + "_MSL.gdb\\"
+Tsunami_2_HAT = output_loc + "Tsunami2-Perth\\" + community + "_HAT.gdb\\"
+Tsunami_2_MSL = output_loc + "Tsunami2-Perth\\" + community + "_MSL.gdb\\"
+Tsunami_3_HAT = output_loc + "Tsunami3-1833\\" + community + "_HAT.gdb\\"
+Tsunami_3_MSL = output_loc + "Tsunami3-1833\\" + community + "_MSL.gdb\\"
+
+hr = community + '_runup'
+
+for location in locations:
+
+    # Local variables...
+    if location == 'Bunbury':
+        # polygon clip
+        extent = input_scenario + "map_work\\Busselton.gdb\\internal_polygons\\bunbury_1km"
+        # name of rasters
+        depth = 'Bun_depth_M_E'
+        speed = 'Bun_speed_M'
+    elif location == 'Busselton':
+        # polygon clip
+        extent = input_scenario + "map_work\\Busselton.gdb\\Internal_polygons\\busselton_1km"
+        # name of rasters
+        depth = 'Bus_depth_M_E'
+        speed = 'Bus_speed_M'
+    else:
+        print 'No recognisable location entered'
+
+    print 'location', location
+    print 'extent', extent
+
+    # Output locations
+    
+    hr_1_HAT = Tsunami_1_HAT +"highest_runup"
+    hr_1_MSL = Tsunami_1_MSL +"highest_runup"
+    hr_2_HAT = Tsunami_2_HAT +"highest_runup"
+    hr_2_MSL = Tsunami_2_MSL +"highest_runup"
+    hr_3_HAT = Tsunami_3_HAT +"highest_runup"
+    hr_3_MSL = Tsunami_3_MSL +"highest_runup"
+
+    depth_1_HAT = Tsunami_1_HAT + location + "_depth"
+    depth_1_MSL = Tsunami_1_MSL + location + "_depth"
+    depth_2_HAT = Tsunami_2_HAT + location + "_depth"
+    depth_2_MSL = Tsunami_2_MSL + location + "_depth"
+    depth_3_HAT = Tsunami_3_HAT + location + "_depth"
+    depth_3_MSL = Tsunami_3_MSL + location + "_depth"
+
+    speed_1_HAT = Tsunami_1_HAT + location + "_speed"
+    speed_1_MSL = Tsunami_1_MSL + location + "_speed"
+    speed_2_HAT = Tsunami_2_HAT + location + "_speed"
+    speed_2_MSL = Tsunami_2_MSL + location + "_speed"
+    speed_3_HAT = Tsunami_3_HAT + location + "_speed"
+    speed_3_MSL = Tsunami_3_MSL + location + "_speed"
+
+
+    # Input locations
+
+    in_depth_1_HAT = in_Tsunami_1_HAT + depth
+    in_depth_1_MSL = in_Tsunami_1_MSL + depth
+    in_depth_2_HAT = in_Tsunami_2_HAT + depth
+    in_depth_2_MSL = in_Tsunami_2_MSL + depth
+    in_depth_3_HAT = in_Tsunami_3_HAT + depth
+    in_depth_3_MSL = in_Tsunami_3_MSL + depth
+
+    in_speed_1_HAT = in_Tsunami_1_HAT + speed
+    in_speed_1_MSL = in_Tsunami_1_MSL + speed
+    in_speed_2_HAT = in_Tsunami_2_HAT + speed
+    in_speed_2_MSL = in_Tsunami_2_MSL + speed
+    in_speed_3_HAT = in_Tsunami_3_HAT + speed
+    in_speed_3_MSL = in_Tsunami_3_MSL + speed
+
+    in_hr_1_HAT = in_Tsunami_1_HAT + hr
+    in_hr_1_MSL = in_Tsunami_1_MSL + hr
+    in_hr_2_HAT = in_Tsunami_2_HAT + hr
+    in_hr_2_MSL = in_Tsunami_2_MSL + hr
+    in_hr_3_HAT = in_Tsunami_3_HAT + hr
+    in_hr_3_MSL = in_Tsunami_3_MSL + hr
+
+
+# Process: Extract by Mask
+    print 'Extract by mask for depth'
+    gp.ExtractByMask_sa(in_depth_1_HAT, extent, depth_1_HAT)
+    gp.ExtractByMask_sa(in_depth_1_MSL, extent, depth_1_MSL)
+    gp.ExtractByMask_sa(in_depth_2_HAT, extent, depth_2_HAT)
+    gp.ExtractByMask_sa(in_depth_2_MSL, extent, depth_2_MSL)
+    gp.ExtractByMask_sa(in_depth_3_HAT, extent, depth_3_HAT)
+    gp.ExtractByMask_sa(in_depth_3_MSL, extent, depth_3_MSL)
+
+    print 'Extract by mask for speed'
+    gp.ExtractByMask_sa(in_speed_1_HAT, extent, speed_1_HAT)
+    gp.ExtractByMask_sa(in_speed_1_MSL, extent, speed_1_MSL)
+    gp.ExtractByMask_sa(in_speed_2_HAT, extent, speed_2_HAT)
+    gp.ExtractByMask_sa(in_speed_2_MSL, extent, speed_2_MSL)
+    gp.ExtractByMask_sa(in_speed_3_HAT, extent, speed_3_HAT)
+    gp.ExtractByMask_sa(in_speed_3_MSL, extent, speed_3_MSL)
+
+
+# Process: Copy Features...
+
+print 'Copy Highest Runup'
+gp.CopyFeatures_management(in_hr_1_HAT, hr_1_HAT, "", "0", "0", "0")
+gp.CopyFeatures_management(in_hr_1_MSL, hr_1_MSL, "", "0", "0", "0")
+gp.CopyFeatures_management(in_hr_2_HAT, hr_2_HAT, "", "0", "0", "0")
+gp.CopyFeatures_management(in_hr_2_MSL, hr_2_MSL, "", "0", "0", "0")
+gp.CopyFeatures_management(in_hr_3_HAT, hr_3_HAT, "", "0", "0", "0")
+gp.CopyFeatures_management(in_hr_3_MSL, hr_3_MSL, "", "0", "0", "0")
+
Index: /anuga_work/production/busselton/busselton_rerun/Inundation_polygon.py
===================================================================
--- /anuga_work/production/busselton/busselton_rerun/Inundation_polygon.py	(revision 7119)
+++ /anuga_work/production/busselton/busselton_rerun/Inundation_polygon.py	(revision 7119)
@@ -0,0 +1,88 @@
+# ---------------------------------------------------------------------------
+# Inundation_line.py
+# Created on: Tue May 12 2009 01:50:45 PM
+#   (generated by ArcGIS/ModelBuilder)
+# Usage: Inundation_line <Bun_depth_M_E> <Bus_depth_M_E> <bunbury_1km> <busselton_1km> 
+# ---------------------------------------------------------------------------
+
+# Import system modules
+import sys, string, os, arcgisscripting, glob
+
+# Create the Geoprocessor object
+gp = arcgisscripting.create()
+gp.overwriteoutput = 1
+
+# Set the necessary product code
+gp.SetProduct("ArcInfo")
+
+# Check out any necessary licenses
+gp.CheckOutExtension("spatial")
+
+# Load required toolboxes...
+gp.AddToolbox("C:/Program Files/ArcGIS/ArcToolbox/Toolboxes/Spatial Analyst Tools.tbx")
+gp.AddToolbox("C:/Program Files/ArcGIS/ArcToolbox/Toolboxes/Conversion Tools.tbx")
+gp.AddToolbox("C:/Program Files/ArcGIS/ArcToolbox/Toolboxes/Data Management Tools.tbx")
+gp.AddToolbox("C:/Program Files/ArcGIS/ArcToolbox/Toolboxes/Analysis Tools.tbx")
+
+# Script arguments...
+time_dir1 = '20090327_152717_run_final_0.6_27283_Bt_kvanputt'
+time_dir2 = '20090330_120458_run_final_0_27283_Bt_kvanputt'
+time_dir3 = '20090324_134328_run_final_0.6_27255_Bt_kvanputt'
+time_dir4 = '20090327_160646_run_final_0_27255_Bt_kvanputt'
+time_dir5 = '20090324_140747_run_final_0.6_68693_Bt_kvanputt'
+time_dir6 = '20090327_161520_run_final_0_68693_Bt_kvanputt'
+
+time_dirs = [time_dir1, time_dir2, time_dir3, time_dir4, time_dir5, time_dir6] 
+
+community = 'busselton'
+scenario = "\\\\nas2\\gemd\\georisk_models\\inundation\\data\\western_australia\\" + community +"_tsunami_scenario\\"
+output_folder = scenario + "anuga\\outputs\\"
+locations = ['bunbury', 'busselton']
+
+for time_dir in time_dirs:
+
+    raster = output_folder + time_dir + "\\raster.gdb\\"
+    
+    for location in locations:
+
+        # Local variables...
+        if location == 'bunbury':
+            # polygon clip
+            extent = scenario + "map_work\\Busselton.gdb\\internal_polygons\\bunbury_1km_clip"
+            # name of rasters
+            depth = raster + 'Bun_depth_M_E'
+            
+        elif location == 'busselton':
+            # polygon clip
+            extent = scenario + "map_work\\Busselton.gdb\\Internal_polygons\\busselton_1km_clip"
+            # name of rasters
+            depth = raster + 'Bus_depth_M_E'
+            
+        else:
+            print 'No recognisable location entered'
+
+        print 'location', location
+        print 'extent', extent
+
+        depth_raster = depth
+
+        # Local variables...
+        Output_raster = raster +"Temp_raster_"+ location
+        Output_polygon_features = raster +"Temp_polygon_"+ location
+        Output_Feature_Class = raster +"Temp_line_"+ location
+        Reclassification = "0 0.05 NODATA;0.05 4.5840001106262207 1"
+        Inundation = raster + "Inundation_" + location + "_polygon"
+        Inundation_line = raster + "Inundation_line"
+          
+        print 'Process: Reclassify'
+        gp.Reclassify_sa(depth, "Value", Reclassification, Output_raster, "DATA")
+        
+        print 'Process: Raster to Polygon'
+        gp.RasterToPolygon_conversion(Output_raster, Output_polygon_features, "SIMPLIFY", "VALUE")
+
+        print 'Process: Clip'
+        gp.Clip_analysis(Output_polygon_features, extent, Inundation, "")
+
+        print 'Process: Delete'
+        gp.Delete_management(Output_raster, "RasterDataset")
+        gp.Delete_management(Output_polygon_features, "FeatureClass")
Index: /anuga_work/production/busselton/busselton_rerun/export_results_max.py
===================================================================
--- /anuga_work/production/busselton/busselton_rerun/export_results_max.py	(revision 7118)
+++ /anuga_work/production/busselton/busselton_rerun/export_results_max.py	(revision 7119)
@@ -30,5 +30,5 @@
 time_dir6 = '20090330_120458_run_final_0_27283_Bt_kvanputt'
 
-time_dirs = [time_dir6] #1,time_dir2, time_dir3, time_dir4, time_dir5]
+time_dirs = [time_dir1] #1,time_dir2, time_dir3, time_dir4, time_dir5]
   
 
@@ -57,5 +57,5 @@
 
 # one or more key strings from var_equations above
-var = ['depth', 'speed']
+var = ['elevation']
 
 ######
@@ -90,7 +90,7 @@
             name3 = join(directory, time_dir, project.scenario_name)+'_time_71040_0'
 
-            names = [name1, name2, name3]
+            names = [name1] #, name2, name3]
       
-            asc_name = []   
+            #asc_name = []   
                
             for name in names:
@@ -113,10 +113,10 @@
                             format = 'asc')
 
-                asc_name.append(outname + '.asc')
-    
-            maxasc_outname = join(directory, time_dir, project.scenario_name)+'_'+which_area+'_'+which_var+'_max.asc'
-
-            print 'max asc outname ', maxasc_outname
-            print 'asc_name ', str(asc_name)
-            
-            MaxAsc(maxasc_outname, asc_name)
+##                asc_name.append(outname + '.asc')
+##    
+##            maxasc_outname = join(directory, time_dir, project.scenario_name)+'_'+which_area+'_'+which_var+'_max.asc'
+##
+##            print 'max asc outname ', maxasc_outname
+##            print 'asc_name ', str(asc_name)
+##            
+##            MaxAsc(maxasc_outname, asc_name)
Index: /anuga_work/production/busselton/busselton_rerun/get_runup_v2.py
===================================================================
--- /anuga_work/production/busselton/busselton_rerun/get_runup_v2.py	(revision 7119)
+++ /anuga_work/production/busselton/busselton_rerun/get_runup_v2.py	(revision 7119)
@@ -0,0 +1,64 @@
+# Standard modules
+from os import sep
+import os
+from os.path import dirname, basename
+from os import mkdir, access, F_OK
+from os.path import join
+from anuga.utilities.polygon import read_polygon
+from shutil import copy
+import time
+import sys
+# Related major packages
+from anuga.shallow_water.data_manager import get_maximum_inundation_data
+from anuga.shallow_water.data_manager import csv2polygons
+
+# Application specific imports
+from setup_model import project  # Definition of file names and polygons
+import csv
+directory = project.output_folder
+
+time_dir1 = '20090324_134328_run_final_0.6_27255_Bt_kvanputt'
+time_dir2 = '20090324_140747_run_final_0.6_68693_Bt_kvanputt'
+time_dir3 = '20090327_152717_run_final_0.6_27283_Bt_kvanputt'
+time_dir4 = '20090327_160646_run_final_0_27255_Bt_kvanputt'
+time_dir5 = '20090327_161520_run_final_0_68693_Bt_kvanputt'
+time_dir6 = '20090330_120458_run_final_0_27283_Bt_kvanputt'
+
+time_dirs = [time_dir1, time_dir4]
+
+
+# polygons is a dictionary of polygon data indexed by 'id'
+# values is a dictionary of polygon names indexed by 'id'
+polygons, values = csv2polygons(project.images)
+print 'Creating %d highest runup' % len(polygons)
+
+# create a dictionary {<polygon_name>: <polygon_data>, ...}
+area_dictionary = {}
+for id in polygons:
+    area_dictionary[values[id]] = polygons[id] 
+
+# for each time directory ...
+for time_dir in time_dirs:
+    print 'time_dir ',time_dir
+    basename = join(directory, time_dir, project.scenario_name)
+
+    out_file = basename+'_runup.csv'
+    output_file = open(out_file, 'w')
+    s = 'area,q,x,y\n' 
+    output_file.write(s)
+
+    # iterate over polygon names
+    for which_area in values.itervalues():
+        print 'which_area = %s' % which_area
+        polygon = area_dictionary[which_area]
+        
+        q, loc = get_maximum_inundation_data(basename+'.sww',
+                                             polygon=polygon, verbose=True)
+
+        print '%s, %.6f, %.6f, %.6f\n' %(str(which_area), q, loc[0],loc[1])
+        s = '%s, %.6f, %.6f, %.6f\n' %(str(which_area), q, loc[0],loc[1])
+
+        output_file.write(s)
+
+        print 'output file printed to: %s' %out_file
+    output_file.close()
Index: /anuga_work/production/busselton/busselton_rerun/project.py
===================================================================
--- /anuga_work/production/busselton/busselton_rerun/project.py	(revision 7118)
+++ /anuga_work/production/busselton/busselton_rerun/project.py	(revision 7119)
@@ -110,4 +110,9 @@
 # Format latitude,longitude etc (geographic)
 building_exposure_filename = 'busselton_res_clip.csv' # from NEXIS
+
+# AREA OF IMAGES - Extent of each image to find out highest runup
+# Header - easting,northing,id,value
+# Used in get_runup.py
+images_filename = 'define_area_clip.csv'
 
 # BOUNDING POLYGON - used in build_boundary.py and run_model.py respectively
@@ -248,4 +253,9 @@
     building_exposure = join(gauges_folder, building_exposure_filename)
 
+# The absolute pathname for the image file
+# Used for get_runup.py
+if images_filename:
+    images = join(polygons_folder, images_filename)
+
 # full path to where MUX files (or meta-files) live
 mux_input = join(event_folder, mux_input_filename)
