Changeset 7568


Ignore:
Timestamp:
Nov 24, 2009, 4:25:22 PM (9 years ago)
Author:
jgriffin
Message:

added code to calculate inundated areas

File:
1 edited

Legend:

Unmodified
Added
Removed
  • anuga_work/production/new_south_wales/batemans_bay/compare_inundation_areas.py

    r7565 r7568  
    99from os.path import join
    1010# Create the Geoprocessor object
    11 gp = arcgisscripting.create()
     11gp = arcgisscripting.create(9.3)
    1212
    1313# Check out any necessary licenses
     
    3737##              'Puysegur_1000yr':58226,
    3838##              'Puysegur_5000yr':58286,
     39##              'Puysegur_15000yr':58272,
    3940##              'New_Hebrides_200yr':51077,
    4041##              'New_Hebrides_500yr':51378,
     
    5354    instantaneous_energy = line_split[2]
    5455    run_up_madsen = line_split[3][0:5]
    55 
     56    event_name = False
    5657    for key, value in event_dict.iteritems():
    5758        if value == int(event_id):
    5859            event_name = key
     60    if event_name == False:
     61        continue
    5962    raster_path = join(filePath,event_name,'raster.gdb')
    60    
     63
    6164   
    6265    print event_id
     
    135138    ####################################################################################
    136139    # Local variables...
    137     inundation = raster_path
     140    gp.workspace = raster_path
     141    print raster_path
     142    inundation_name = gp.ListRasters("*depth*","")
     143    for raster in inundation_name:
     144        print raster
     145    inundation = raster_path + "\\"+ inundation_name[0]
    138146    inundation_area1 = "N:\\georisk_models\\inundation\data\\new_south_wales\\batemans_bay_tsunami_scenario_2009\\anuga\\polygons\\polygons.gdb\inundation_area1"
    139     inundation_reclass = filePath+"elevation\\elevation.gdb\\elevation_"+run_up_name+"_reclass"
    140     inundation_poly = filePath+"elevation\\features.gdb\\elevation_"+run_up_name+"_poly"
    141     inundation_less = filePath+"elevation\\features.gdb\\elevation_"+run_up_name+"_less"
    142     inundation_less_final = filePath+"elevation\\features.gdb\\elevation_less_final"
    143 
    144     # Process: Extract elevation data by Mask...
    145     print 'check if elevation already extracted'
    146     if gp.Exists(elevation_extract)!=True:
     147    inundation_extract = "inundation_"+run_up_name+"_extract"
     148    inundation_reclass = "inundation_"+run_up_name+"_reclass"
     149    inundation_poly = "inundation_"+run_up_name+"_poly"
     150    inundation_less = "inundation_"+run_up_name+"_less"
     151    inundation_less_final = "inundation_less_final"
     152
     153    # Process: Extract inundation data by Mask...
     154    print 'check if inundation already extracted'
     155    if gp.Exists(inundation_extract)!=True:
    147156        print 'not already extracted'
    148157        print 'extracting by mask'
    149         gp.ExtractByMask_sa(elevation, inundation_area1, elevation_extract)
    150     else:
    151         print 'extracted elevation already exists'
    152 
    153     # Process: Reclassify elevation data
    154     print 'check if elevation already reclassified'
    155     if gp.Exists(elevation_reclass)!=True:
    156         print 'reclassify based on elevation height of:', run_up_height
    157         reclass_level = "-200 "+run_up_height+" 1;"+run_up_height+" 300 0"
    158         gp.Reclassify_sa(elevation_extract, "Value", reclass_level, elevation_reclass, "DATA")
    159     else:
    160         print 'elevation has previously been reclassified for this height:', run_up_height
     158        gp.ExtractByMask_sa(inundation, inundation_area1, inundation_extract)
     159    else:
     160        print 'extracted inundation already exists'
     161
     162    # Process: Reclassify inundation data
     163    print 'check if inundation already reclassified'
     164    if gp.Exists(inundation_reclass)!=True:
     165        print 'reclassify based on inundation'
     166        gp.Reclassify_sa(inundation_extract, "Value", "-200 0 1;0 300 0", inundation_reclass, "DATA")
     167    else:
     168        print 'inundation has previously been reclassified for this height:', run_up_height
    161169
    162170    # Process: Raster to Polygon...
    163171    print 'check if already converted from raster to polyon'
    164     if gp.Exists(elevation_poly)!=True:
     172    if gp.Exists(inundation_poly)!=True:
    165173        print 'raster to polyon'
    166         gp.RasterToPolygon_conversion(elevation_reclass, elevation_poly, "NO_SIMPLIFY", "VALUE")
    167     else:
    168         print 'elevation raster already converted to polygon'
     174        gp.RasterToPolygon_conversion(inundation_reclass, inundation_poly, "NO_SIMPLIFY", "VALUE")
     175    else:
     176        print 'inundation raster already converted to polygon'
    169177
    170178    # Process: Select...
    171179    print 'select by attribute'
    172     gp.Select_analysis(elevation_poly, elevation_less, "\"GRIDCODE\" = 1")
     180    gp.Select_analysis(inundation_poly, inundation_less, "\"GRIDCODE\" = 1")
    173181
    174182    # Process: Create layer...
    175183    print 'creating layers'
    176184    gp.MakeFeatureLayer(batemans_bay_SMARTLINE,'smartline_layer')
    177     gp.MakeFeatureLayer(elevation_less,'elevation_layer')
     185    gp.MakeFeatureLayer(inundation_less,'inundation_layer')
    178186
    179187    # Process: Select Layer By Location...
    180188    print 'select by location'
    181     gp.SelectLayerByLocation_management('elevation_layer', "INTERSECT", 'smartline_layer', "", "NEW_SELECTION")
     189    gp.SelectLayerByLocation_management('inundation_layer', "INTERSECT", 'smartline_layer', "", "NEW_SELECTION")
    182190    print 'joining'
    183     gp.SpatialJoin_analysis('elevation_layer',inundation_area1, elevation_less_final,"JOIN_ONE_TO_ONE","","","INTERSECTS")
     191    gp.SpatialJoin_analysis('inundation_layer',inundation_area1, inundation_less_final,"JOIN_ONE_TO_ONE","","","INTERSECTS")
    184192
    185193
Note: See TracChangeset for help on using the changeset viewer.