Changeset 7565


Ignore:
Timestamp:
Nov 23, 2009, 4:08:17 PM (15 years ago)
Author:
jgriffin
Message:
 
File:
1 edited

Legend:

Unmodified
Added
Removed

  • anuga_work/production/new_south_wales/batemans_bay/compare_inundation_areas.py

    r7563 r7565  
    2929event_dict = {'Event1_MSL':58346,
    3030              'Event1_HAT':58346,
    31               'Event2_MSL':51204,
    32               'Event2_HAT':51204,
    33               'Event3_MSL':58284,
    34               'Event3_HAT':58284,
    35               'Puysegur_200yr':58129,
    36               'Puysegur_500yr':58115,
    37               'Puysegur_1000yr':58226,
    38               'Puysegur_5000yr':58286,
    39               'New_Hebrides_200yr':51077,
    40               'New_Hebrides_500yr':51378,
    41               'New_Hebrides_1000yr':51347,
    42               'New_Hebrides_2000yr':51292,
    43               'New_Hebrides_5000yr':51424,}
     31##              'Event2_MSL':51204,
     32##              'Event2_HAT':51204,
     33              'Event3_MSL':58284}#,
     34##              'Event3_HAT':58284,
     35##              'Puysegur_200yr':58129,
     36##              'Puysegur_500yr':58115,
     37##              'Puysegur_1000yr':58226,
     38##              'Puysegur_5000yr':58286,
     39##              'New_Hebrides_200yr':51077,
     40##              'New_Hebrides_500yr':51378,
     41##              'New_Hebrides_1000yr':51347,
     42##              'New_Hebrides_2000yr':51292,
     43##              'New_Hebrides_5000yr':51424,}
    4444
    4545boundaries = file(boundary_file,'r')
    … …  
    5050    gauge_file = line_split[0]
    5151    event_id = gauge_file.split('/')[10]
     52    crest_integrated_energy = line_split[1]
     53    instantaneous_energy = line_split[2]
     54    run_up_madsen = line_split[3][0:5]
     55
     56    for key, value in event_dict.iteritems():
     57        if value == int(event_id):
     58            event_name = key
     59    raster_path = join(filePath,event_name,'raster.gdb')
    5260   
    53     break
    5461   
    55 
    56 
    57 
    58 
    59 
     62    print event_id
     63    print gauge_file
     64    #break
     65
     66
     67
     68   
     69    ##run_up_madsen = 6.0
     70    run_up_height = str(run_up_madsen)
     71    run_up_name = run_up_height.split('.')[0]
     72
     73
     74    ####################################################################################
     75    ##    Extract areas for analytical solution
     76    ####################################################################################
    6077    # Local variables...
    61     elevation_6m_split = filePath+"elevation\\elevation.gdb\\elevation_6m_split"
    62     elevation_6m_split__2_ = filePath+"elevation\\features.gdb\\elevation_6m_split"
    6378    elevation = filePath+"elevation\\elevation.gdb\\elevation"
    64     elevation_6m_less = filePath+"elevation\\features.gdb\\elevation_6m_less"
    6579    elevation_extract = filePath+"elevation\\elevation.gdb\\elevation_extract"
    6680    inundation_area1 = "N:\\georisk_models\\inundation\data\\new_south_wales\\batemans_bay_tsunami_scenario_2009\\anuga\\polygons\\polygons.gdb\inundation_area1"
    67     elevation_6m_less__2_ = filePath+"elevation\\features.gdb\\elevation_6m_less"
    68     elevation_6m_less__3_ = filePath+"elevation\\features.gdb\\elevation_6m_less"
     81    elevation_reclass = filePath+"elevation\\elevation.gdb\\elevation_"+run_up_name+"_reclass"
     82    elevation_poly = filePath+"elevation\\features.gdb\\elevation_"+run_up_name+"_poly"
     83    elevation_less = filePath+"elevation\\features.gdb\\elevation_"+run_up_name+"_less"
    6984    batemans_bay_SMARTLINE = "N:\\georisk_models\\inundation\\data\\new_south_wales\\coastline\\batemans_bay_SMARTLINE.shp"
    70     elevation_6m_less_final = filePath+"elevation\\features.gdb\\elevation_6m_less_final"
    71 
    72     # Process: Extract by Mask...
    73     print 'extracting by mask'
    74     gp.ExtractByMask_sa(elevation, inundation_area1, elevation_extract)
    75 
    76     # Process: Reclassify...
    77     print 'reclassify based on elevation height'
    78     gp.Reclassify_sa(elevation_extract, "Value", "-200 6 1;6 300 0", elevation_6m_split, "DATA")
     85    elevation_less_final = filePath+"elevation\\features.gdb\\elevation_less_final"
     86
     87    # Process: Extract elevation data by Mask...
     88    print 'check if elevation already extracted'
     89    if gp.Exists(elevation_extract)!=True:
     90        print 'not already extracted'
     91        print 'extracting by mask'
     92        gp.ExtractByMask_sa(elevation, inundation_area1, elevation_extract)
     93    else:
     94        print 'extracted elevation already exists'
     95
     96    # Process: Reclassify elevation data
     97    print 'check if elevation already reclassified'
     98    if gp.Exists(elevation_reclass)!=True:
     99        print 'reclassify based on elevation height of:', run_up_height
     100        reclass_level = "-200 "+run_up_height+" 1;"+run_up_height+" 300 0"
     101        gp.Reclassify_sa(elevation_extract, "Value", reclass_level, elevation_reclass, "DATA")
     102    else:
     103        print 'elevation has previously been reclassified for this height:', run_up_height
    79104
    80105    # Process: Raster to Polygon...
    81     print 'raster to polyon'
    82     gp.RasterToPolygon_conversion(elevation_6m_split, elevation_6m_split__2_, "NO_SIMPLIFY", "VALUE")
     106    print 'check if already converted from raster to polyon'
     107    if gp.Exists(elevation_poly)!=True:
     108        print 'raster to polyon'
     109        gp.RasterToPolygon_conversion(elevation_reclass, elevation_poly, "NO_SIMPLIFY", "VALUE")
     110    else:
     111        print 'elevation raster already converted to polygon'
    83112
    84113    # Process: Select...
    85114    print 'select by attribute'
    86     gp.Select_analysis(elevation_6m_split__2_, elevation_6m_less, "\"GRIDCODE\" = 1")
     115    gp.Select_analysis(elevation_poly, elevation_less, "\"GRIDCODE\" = 1")
    87116
    88117    # Process: Create layer...
    89     print 'creating layer'
     118    print 'creating layers'
    90119    gp.MakeFeatureLayer(batemans_bay_SMARTLINE,'smartline_layer')
    91     gp.MakeFeatureLayer(elevation_6m_less__3_,'elevation_layer')
     120    gp.MakeFeatureLayer(elevation_less,'elevation_layer')
    92121
    93122    # Process: Select Layer By Location...
    … …  
    95124    gp.SelectLayerByLocation_management('elevation_layer', "INTERSECT", 'smartline_layer', "", "NEW_SELECTION")
    96125    print 'joining'
    97     gp.SpatialJoin_analysis('elevation_layer',inundation_area1, elevation_6m_less_final,"JOIN_ONE_TO_ONE","","","INTERSECTS")
     126    gp.SpatialJoin_analysis('elevation_layer',inundation_area1, elevation_less_final,"JOIN_ONE_TO_ONE","","","INTERSECTS")
    98127
    99128    # Process: Copy Features...
    100129    ##print 'copy features to output feature class'
    101     ##gp.CopyFeatures_management('elevation_layer', elevation_6m_less_final, "", "0", "0", "0")
     130    ##gp.CopyFeatures_management('elevation_layer', elevation_less_final, "", "0", "0", "0")
     131
     132
     133    ####################################################################################
     134    ##    Extract areas for ANUGA solution
     135    ####################################################################################
     136    # Local variables...
     137    inundation = raster_path
     138    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        print 'not already extracted'
     148        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
     161
     162    # Process: Raster to Polygon...
     163    print 'check if already converted from raster to polyon'
     164    if gp.Exists(elevation_poly)!=True:
     165        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'
     169
     170    # Process: Select...
     171    print 'select by attribute'
     172    gp.Select_analysis(elevation_poly, elevation_less, "\"GRIDCODE\" = 1")
     173
     174    # Process: Create layer...
     175    print 'creating layers'
     176    gp.MakeFeatureLayer(batemans_bay_SMARTLINE,'smartline_layer')
     177    gp.MakeFeatureLayer(elevation_less,'elevation_layer')
     178
     179    # Process: Select Layer By Location...
     180    print 'select by location'
     181    gp.SelectLayerByLocation_management('elevation_layer', "INTERSECT", 'smartline_layer', "", "NEW_SELECTION")
     182    print 'joining'
     183    gp.SpatialJoin_analysis('elevation_layer',inundation_area1, elevation_less_final,"JOIN_ONE_TO_ONE","","","INTERSECTS")
     184
     185
    102186
    103187
    … …  
    105189    # Search feature class for areas
    106190    print 'get areas'
     191    area_dict = {}
    107192    areas = []
    108     cur = gp.SearchCursor(elevation_6m_less_final)
     193    cur = gp.SearchCursor(elevation_less_final)
    109194    sRow = cur.Next()
    110195    while sRow:
     196        area_dict[sRow.GetValue("inundation_ID")] = sRow.GetValue("Shape_Area")
    111197        areas.append(sRow.GetValue("Shape_Area"))
    112198        sRow = cur.Next()
    113199    print areas
     200    print area_dict
     201
     202boundaries.close()
Note: See TracChangeset for help on using the changeset viewer.