Changeset 7570

Timestamp:

Nov 25, 2009, 4:28:13 PM (13 years ago)

Author:

jgriffin

Message:

 

File:

• anuga_work/production/new_south_wales/batemans_bay/compare_inundation_areas.py (modified) (5 diffs)

anuga_work/production/new_south_wales/batemans_bay/compare_inundation_areas.py

@@ +1 @@
+"""# ---------------------------------------------------------------------------
+ elevation_model.py
+ Created on: Tue Nov 17 2009 02:11:08 
+ Created by: Jonathan Griffin
+ This script compares calculated inundate areas for run-up solutions from analystical models
+ and ANUGA modelling.
+
+
 # ---------------------------------------------------------------------------
-# 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
+from collections import defaultdict
 # Create the Geoprocessor object
 gp = arcgisscripting.create(9.3)
@@ -43 +48 @@
 ##              'New_Hebrides_2000yr':51292,
 ##              'New_Hebrides_5000yr':51424,}
+
+# Dictionary for storing area comparison
+area_comp = defaultdict(list)
 
 boundaries = file(boundary_file,'r')
@@ -72 +80 @@
     ##run_up_madsen = 6.0
     run_up_height = str(run_up_madsen)
-    run_up_name = run_up_height.split('.')[0]
+    run_up_name = run_up_height.replace('.','_')[0]
 
 
@@ -178 +186 @@
     # Process: Select...
     print 'select by attribute'
-    gp.Select_analysis(inundation_poly, inundation_less, "\"GRIDCODE\" = 1")
+    gp.Select_analysis(inundation_poly, inundation_less, "\"GRIDCODE\" = 0")
 
     # Process: Create layer...
@@ -197 +205 @@
     # Search feature class for areas
     print 'get areas'
-    area_dict = {}
-    areas = []
+    analytical_area_dict = {}
+    #analytical_areas = []
     cur = gp.SearchCursor(elevation_less_final)
     sRow = cur.Next()
     while sRow:
-        area_dict[sRow.GetValue("inundation_ID")] = sRow.GetValue("Shape_Area")
-        areas.append(sRow.GetValue("Shape_Area"))
+        if not(sRow.GetValue("inundation_ID")in analytical_area_dict):
+            analytical_area_dict[sRow.GetValue("inundation_ID")] = sRow.GetValue("Shape_Area")
+        else:
+            analytical_area_dict[sRow.GetValue("inundation_ID")] = (analytical_area_dict[sRow.GetValue("inundation_ID")]+sRow.GetValue("Shape_Area"))
+        #analytical_areas.append(sRow.GetValue("Shape_Area"))
         sRow = cur.Next()
-    print areas
-    print area_dict
-
+    #print analytical_areas
+    print analytical_area_dict
+
+    # Search feature class for areas
+    print 'get areas'
+    anuga_area_dict = {}
+    #anuga_areas = []
+    cur = gp.SearchCursor(inundation_less_final)
+    sRow = cur.Next()
+    while sRow:
+        if not(sRow.GetValue("inundation_ID") in anuga_area_dict):
+            anuga_area_dict[sRow.GetValue("inundation_ID")]=sRow.GetValue("Shape_Area")
+        else:
+            anuga_area_dict[sRow.GetValue("inundation_ID")]= anuga_area_dict[sRow.GetValue("inundation_ID")]+sRow.GetValue("Shape_Area")  
+        #anuga_areas.append(sRow.GetValue("Shape_Area"))
+        sRow = cur.Next()
+    #print anuga_areas
+    print anuga_area_dict
+
+    for key in anuga_area_dict:
+        diff = anuga_area_dict[key]-analytical_area_dict[key]
+        area_comp[key].append(diff)
+    print 'area_comp', area_comp
+    
 boundaries.close()