Changeset 7568

Timestamp:

Nov 24, 2009, 4:25:22 PM (15 years ago)

Author:

jgriffin

Message:

added code to calculate inundated areas

File:

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

anuga_work/production/new_south_wales/batemans_bay/compare_inundation_areas.py

@@ -9 +9 @@
 from os.path import join
 # Create the Geoprocessor object
-gp = arcgisscripting.create()
+gp = arcgisscripting.create(9.3)
 
 # Check out any necessary licenses
@@ -37 +37 @@
 ##              'Puysegur_1000yr':58226,
 ##              'Puysegur_5000yr':58286,
+##              'Puysegur_15000yr':58272,
 ##              'New_Hebrides_200yr':51077,
 ##              'New_Hebrides_500yr':51378,
@@ -53 +54 @@
     instantaneous_energy = line_split[2]
     run_up_madsen = line_split[3][0:5]
-
+    event_name = False
     for key, value in event_dict.iteritems():
         if value == int(event_id):
             event_name = key
+    if event_name == False:
+        continue
     raster_path = join(filePath,event_name,'raster.gdb')
-    
+
     
     print event_id
@@ -135 +138 @@
     ####################################################################################
     # Local variables...
-    inundation = raster_path
+    gp.workspace = raster_path
+    print raster_path
+    inundation_name = gp.ListRasters("*depth*","")
+    for raster in inundation_name:
+        print raster
+    inundation = raster_path + "\\"+ inundation_name[0]
     inundation_area1 = "N:\\georisk_models\\inundation\data\\new_south_wales\\batemans_bay_tsunami_scenario_2009\\anuga\\polygons\\polygons.gdb\inundation_area1"
-    inundation_reclass = filePath+"elevation\\elevation.gdb\\elevation_"+run_up_name+"_reclass"
-    inundation_poly = filePath+"elevation\\features.gdb\\elevation_"+run_up_name+"_poly"
-    inundation_less = filePath+"elevation\\features.gdb\\elevation_"+run_up_name+"_less"
-    inundation_less_final = filePath+"elevation\\features.gdb\\elevation_less_final"
-
-    # Process: Extract elevation data by Mask...
-    print 'check if elevation already extracted'
-    if gp.Exists(elevation_extract)!=True:
+    inundation_extract = "inundation_"+run_up_name+"_extract"
+    inundation_reclass = "inundation_"+run_up_name+"_reclass"
+    inundation_poly = "inundation_"+run_up_name+"_poly"
+    inundation_less = "inundation_"+run_up_name+"_less"
+    inundation_less_final = "inundation_less_final"
+
+    # Process: Extract inundation data by Mask...
+    print 'check if inundation already extracted'
+    if gp.Exists(inundation_extract)!=True:
         print 'not already extracted'
         print 'extracting by mask'
-        gp.ExtractByMask_sa(elevation, inundation_area1, elevation_extract)
-    else:
-        print 'extracted elevation already exists'
-
-    # Process: Reclassify elevation data
-    print 'check if elevation already reclassified'
-    if gp.Exists(elevation_reclass)!=True:
-        print 'reclassify based on elevation height of:', run_up_height
-        reclass_level = "-200 "+run_up_height+" 1;"+run_up_height+" 300 0"
-        gp.Reclassify_sa(elevation_extract, "Value", reclass_level, elevation_reclass, "DATA")
-    else:
-        print 'elevation has previously been reclassified for this height:', run_up_height
+        gp.ExtractByMask_sa(inundation, inundation_area1, inundation_extract)
+    else:
+        print 'extracted inundation already exists'
+
+    # Process: Reclassify inundation data
+    print 'check if inundation already reclassified'
+    if gp.Exists(inundation_reclass)!=True:
+        print 'reclassify based on inundation'
+        gp.Reclassify_sa(inundation_extract, "Value", "-200 0 1;0 300 0", inundation_reclass, "DATA")
+    else:
+        print 'inundation has previously been reclassified for this height:', run_up_height
 
     # Process: Raster to Polygon...
     print 'check if already converted from raster to polyon'
-    if gp.Exists(elevation_poly)!=True:
+    if gp.Exists(inundation_poly)!=True:
         print 'raster to polyon'
-        gp.RasterToPolygon_conversion(elevation_reclass, elevation_poly, "NO_SIMPLIFY", "VALUE")
-    else:
-        print 'elevation raster already converted to polygon'
+        gp.RasterToPolygon_conversion(inundation_reclass, inundation_poly, "NO_SIMPLIFY", "VALUE")
+    else:
+        print 'inundation raster already converted to polygon'
 
     # Process: Select...
     print 'select by attribute'
-    gp.Select_analysis(elevation_poly, elevation_less, "\"GRIDCODE\" = 1")
+    gp.Select_analysis(inundation_poly, inundation_less, "\"GRIDCODE\" = 1")
 
     # Process: Create layer...
     print 'creating layers'
     gp.MakeFeatureLayer(batemans_bay_SMARTLINE,'smartline_layer')
-    gp.MakeFeatureLayer(elevation_less,'elevation_layer')
+    gp.MakeFeatureLayer(inundation_less,'inundation_layer')
 
     # Process: Select Layer By Location...
     print 'select by location'
-    gp.SelectLayerByLocation_management('elevation_layer', "INTERSECT", 'smartline_layer', "", "NEW_SELECTION")
+    gp.SelectLayerByLocation_management('inundation_layer', "INTERSECT", 'smartline_layer', "", "NEW_SELECTION")
     print 'joining'
-    gp.SpatialJoin_analysis('elevation_layer',inundation_area1, elevation_less_final,"JOIN_ONE_TO_ONE","","","INTERSECTS")
+    gp.SpatialJoin_analysis('inundation_layer',inundation_area1, inundation_less_final,"JOIN_ONE_TO_ONE","","","INTERSECTS")