Index: /anuga_work/production/new_south_wales/batemans_bay/compare_inundation_areas.py
===================================================================
--- /anuga_work/production/new_south_wales/batemans_bay/compare_inundation_areas.py	(revision 7563)
+++ /anuga_work/production/new_south_wales/batemans_bay/compare_inundation_areas.py	(revision 7563)
@@ -0,0 +1,113 @@
+# ---------------------------------------------------------------------------
+# elevation_model.py
+# Created on: Tue Nov 17 2009 02:11:08 
+#   (generated by ArcGIS/ModelBuilder)
+# ---------------------------------------------------------------------------
+
+# Import system modules
+import sys, string, os, arcgisscripting
+from os.path import join
+# 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/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")
+
+gp.overwriteoutput = 1
+
+filePath = "N:\\georisk_models\\inundation\\data\\new_south_wales\\batemans_bay_tsunami_scenario_2009\\anuga\\outputs\\"
+boundary_path = "N:\\georisk_models\\inundation\\data\\new_south_wales\\batemans_bay_tsunami_scenario_2009\\anuga\\boundaries\\"
+boundary_file_name = 'wave_energy_summary.csv'
+boundary_file = join(boundary_path,boundary_file_name)
+
+event_dict = {'Event1_MSL':58346,
+              'Event1_HAT':58346,
+              'Event2_MSL':51204,
+              'Event2_HAT':51204,
+              'Event3_MSL':58284,
+              'Event3_HAT':58284,
+              'Puysegur_200yr':58129,
+              'Puysegur_500yr':58115,
+              'Puysegur_1000yr':58226,
+              'Puysegur_5000yr':58286,
+              'New_Hebrides_200yr':51077,
+              'New_Hebrides_500yr':51378,
+              'New_Hebrides_1000yr':51347,
+              'New_Hebrides_2000yr':51292,
+              'New_Hebrides_5000yr':51424,}
+
+boundaries = file(boundary_file,'r')
+for line in boundaries.readlines():
+    line_split = line.split(',')
+    if line_split[0]=='filename':
+        continue
+    gauge_file = line_split[0]
+    event_id = gauge_file.split('/')[10]
+    
+    break
+    
+
+
+
+
+
+    # Local variables...
+    elevation_6m_split = filePath+"elevation\\elevation.gdb\\elevation_6m_split"
+    elevation_6m_split__2_ = filePath+"elevation\\features.gdb\\elevation_6m_split"
+    elevation = filePath+"elevation\\elevation.gdb\\elevation"
+    elevation_6m_less = filePath+"elevation\\features.gdb\\elevation_6m_less"
+    elevation_extract = filePath+"elevation\\elevation.gdb\\elevation_extract"
+    inundation_area1 = "N:\\georisk_models\\inundation\data\\new_south_wales\\batemans_bay_tsunami_scenario_2009\\anuga\\polygons\\polygons.gdb\inundation_area1"
+    elevation_6m_less__2_ = filePath+"elevation\\features.gdb\\elevation_6m_less"
+    elevation_6m_less__3_ = filePath+"elevation\\features.gdb\\elevation_6m_less"
+    batemans_bay_SMARTLINE = "N:\\georisk_models\\inundation\\data\\new_south_wales\\coastline\\batemans_bay_SMARTLINE.shp"
+    elevation_6m_less_final = filePath+"elevation\\features.gdb\\elevation_6m_less_final"
+
+    # Process: Extract by Mask...
+    print 'extracting by mask'
+    gp.ExtractByMask_sa(elevation, inundation_area1, elevation_extract)
+
+    # Process: Reclassify...
+    print 'reclassify based on elevation height'
+    gp.Reclassify_sa(elevation_extract, "Value", "-200 6 1;6 300 0", elevation_6m_split, "DATA")
+
+    # Process: Raster to Polygon...
+    print 'raster to polyon'
+    gp.RasterToPolygon_conversion(elevation_6m_split, elevation_6m_split__2_, "NO_SIMPLIFY", "VALUE")
+
+    # Process: Select...
+    print 'select by attribute'
+    gp.Select_analysis(elevation_6m_split__2_, elevation_6m_less, "\"GRIDCODE\" = 1")
+
+    # Process: Create layer...
+    print 'creating layer'
+    gp.MakeFeatureLayer(batemans_bay_SMARTLINE,'smartline_layer')
+    gp.MakeFeatureLayer(elevation_6m_less__3_,'elevation_layer')
+
+    # Process: Select Layer By Location...
+    print 'select by location'
+    gp.SelectLayerByLocation_management('elevation_layer', "INTERSECT", 'smartline_layer', "", "NEW_SELECTION")
+    print 'joining'
+    gp.SpatialJoin_analysis('elevation_layer',inundation_area1, elevation_6m_less_final,"JOIN_ONE_TO_ONE","","","INTERSECTS")
+
+    # Process: Copy Features...
+    ##print 'copy features to output feature class'
+    ##gp.CopyFeatures_management('elevation_layer', elevation_6m_less_final, "", "0", "0", "0")
+
+
+
+    # Search feature class for areas
+    print 'get areas'
+    areas = []
+    cur = gp.SearchCursor(elevation_6m_less_final)
+    sRow = cur.Next()
+    while sRow:
+        areas.append(sRow.GetValue("Shape_Area"))
+        sRow = cur.Next()
+    print areas
