Changeset 7592


Ignore:
Timestamp:
Dec 17, 2009, 9:35:51 AM (14 years ago)
Author:
jgriffin
Message:
 
Location:
anuga_work/production/new_south_wales/batemans_bay
Files:
1 added
5 edited

Legend:

Unmodified
Added
Removed

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

    r7577 r7592  
    3535boundary_file_name = 'wave_energy_summary.csv'
    3636out_file_name_madsen = 'area_comparisons_madsen.csv'
     37out_file_name_surf_sim = 'area_comparisons_surf_sim.csv'
    3738out_file_name_energy = 'area_comparisons_energy.csv'
    3839out_file_name_anuga = 'area_comparisons_anuga.csv'
    … …  
    4243boundary_file = join(boundary_path,boundary_file_name)
    4344out_file_madsen = join(boundary_path, out_file_name_madsen)
     45out_file_surf_sim = join(boundary_path, out_file_name_surf_sim)
    4446out_file_energy = join(boundary_path, out_file_name_energy)
    4547out_file_anuga = join(boundary_path, out_file_name_anuga)
    … …  
    5456##              'Event3_HAT':58284,
    5557              'Puysegur_200yr':58129,
    56 ##              'Puysegur_500yr':58115,
     58              'Puysegur_500yr':58115,
    5759              'Puysegur_1000yr':58226,
    5860              'Puysegur_5000yr':58286,
    … …  
    6466              'New_Hebrides_5000yr':51424,}
    6567
     68# Set True to skip over features already created, False to recreate all features
     69do_all = False
     70
    6671# Dictionaries for storing area comparison
    6772area_comp_madsen = defaultdict(list)
    6873area_comp_energy = defaultdict(list)
    6974area_percent_madsen = defaultdict(list)
     75area_comp_surf_sim = defaultdict(list)
     76area_percent_surf_sim = defaultdict(list)
    7077area_percent_energy = defaultdict(list)
    7178area_total_anuga = defaultdict(list)
    … …  
    8390    run_up_madsen = line_split[3][0:5]
    8491    run_up_energy = line_split[4][0:5]
     92    run_up_surf_sim = line_split[6][0:5]
    8593    event_name = False
    8694    for key, value in event_dict.iteritems():
    … …  
    8997    if event_name == False:
    9098        continue
    91     raster_path = join(filePath,event_name,'raster.gdb')
     99    raster_path = join(filePath,event_name,'area_comparisons.gdb')
    92100
    93101   
    … …  
    96104    #break
    97105
    98     do_all = True
     106   
    99107
    100108   
    … …  
    104112    energy_run_up_height = str(run_up_energy)
    105113    energy_run_up_name = run_up_energy.replace('.','_')[0:5]
     114    surf_sim_name = run_up_surf_sim.replace('.','_')[0:5]
    106115
    107116
    … …  
    124133##    if do_all!=True:
    125134##        print 'continuing'
    126     if gp.Exists(elevation_less_final_dis):
     135    if (gp.Exists(elevation_less_final_dis) and do_all == False):
    127136        print 'analysis already done, skipping to next one'
    128137    else:
    129138        # Process: Extract elevation data by Mask...
    130139        print 'check if elevation already extracted'
    131         if gp.Exists(elevation_extract)!=True:
     140        if (gp.Exists(elevation_extract)!=True or do_all == True):
    132141            print 'not already extracted'
    133142            print 'extracting by mask'
    … …  
    138147        # Process: Reclassify elevation data
    139148        print 'check if elevation already reclassified'
    140         if gp.Exists(elevation_reclass)!=True:
     149        if (gp.Exists(elevation_reclass)!=True or do_all == True):
    141150            print 'reclassify based on elevation height of:', run_up_height
    142151            reclass_level = "-200 "+run_up_height+" 1;"+run_up_height+" 300 0"
    … …  
    147156        # Process: Raster to Polygon...
    148157        print 'check if already converted from raster to polyon'
    149         if gp.Exists(elevation_poly)!=True:
     158        if (gp.Exists(elevation_poly)!=True or do_all == True):
    150159            print 'raster to polyon'
    151160            gp.RasterToPolygon_conversion(elevation_reclass, elevation_poly, "NO_SIMPLIFY", "VALUE")
    … …  
    173182        gp.Dissolve_management(elevation_less_final,elevation_less_final_dis,'inundation_ID')
    174183
    175 
     184    ####################################################################################
     185    ##    Extract areas for Madsen surf similarity analytical solution
     186    ####################################################################################
     187    print '\nDoing Madsen surf_similarity solution...'
     188    # Local variables...
     189    elevation = filePath+"elevation\\elevation.gdb\\elevation"
     190    elevation_extract = filePath+"elevation\\elevation.gdb\\elevation_extract"
     191    inundation_area1 = "N:\\georisk_models\\inundation\data\\new_south_wales\\batemans_bay_tsunami_scenario_2009\\anuga\\polygons\\polygons.gdb\inundation_area1"
     192    elevation_reclass = filePath+"elevation\\elevation.gdb\\elevation_"+surf_sim_name+"_reclass"
     193    elevation_poly = filePath+"elevation\\features.gdb\\elevation_"+surf_sim_name+"_poly"
     194    elevation_less = filePath+"elevation\\features.gdb\\elevation_"+surf_sim_name+"_less"
     195    batemans_bay_SMARTLINE = "N:\\georisk_models\\inundation\\data\\new_south_wales\\coastline\\batemans_bay_SMARTLINE.shp"
     196    elevation_surf_sim_less_final = filePath+"elevation\\features.gdb\\elevation_"+surf_sim_name+"_less_final"
     197    elevation_surf_sim_less_final_dis = filePath+"elevation\\features.gdb\\elevation_"+surf_sim_name+"_less_final_dis"
     198
     199    # Check if already done...
     200##    if do_all!=True:
     201##        print 'continuing'
     202    if (gp.Exists(elevation_surf_sim_less_final_dis) and do_all == False):
     203        print 'analysis already done, skipping to next one'
     204    else:
     205        # Process: Extract elevation data by Mask...
     206        print 'check if elevation already extracted'
     207        if (gp.Exists(elevation_extract)!=True or do_all == True):
     208            print 'not already extracted'
     209            print 'extracting by mask'
     210            gp.ExtractByMask_sa(elevation, inundation_area1, elevation_extract)
     211        else:
     212            print 'extracted elevation already exists'
     213
     214        # Process: Reclassify elevation data
     215        print 'check if elevation already reclassified'
     216        if (gp.Exists(elevation_reclass)!=True or do_all == True):
     217            print 'reclassify based on elevation height of:', run_up_height
     218            reclass_level = "-200 "+run_up_height+" 1;"+run_up_height+" 300 0"
     219            gp.Reclassify_sa(elevation_extract, "Value", reclass_level, elevation_reclass, "DATA")
     220        else:
     221            print 'elevation has previously been reclassified for this height:', run_up_height
     222
     223        # Process: Raster to Polygon...
     224        print 'check if already converted from raster to polyon'
     225        if (gp.Exists(elevation_poly)!=True or do_all == True):
     226            print 'raster to polyon'
     227            gp.RasterToPolygon_conversion(elevation_reclass, elevation_poly, "NO_SIMPLIFY", "VALUE")
     228        else:
     229            print 'elevation raster already converted to polygon'
     230
     231        # Process: Select...
     232        print 'select by attribute'
     233        gp.Select_analysis(elevation_poly, elevation_less, "\"GRIDCODE\" = 1")
     234
     235        # Process: Create layer...
     236        print 'creating layers'
     237        gp.MakeFeatureLayer(batemans_bay_SMARTLINE,'smartline_layer')
     238        gp.MakeFeatureLayer(elevation_less,'elevation_layer')
     239
     240        # Process: Select Layer By Location...
     241        print 'select by location'
     242        gp.SelectLayerByLocation_management('elevation_layer', "INTERSECT", 'smartline_layer', "", "NEW_SELECTION")
     243        print 'joining'
     244        print 'inundation_area1',inundation_area1
     245        print 'energy_less_final',elevation_surf_sim_less_final
     246        gp.SpatialJoin_analysis('elevation_layer',inundation_area1, elevation_surf_sim_less_final,"JOIN_ONE_TO_ONE","","","INTERSECTS")
     247
     248        print 'dissolving fields'
     249        gp.Dissolve_management(elevation_surf_sim_less_final,elevation_surf_sim_less_final_dis,'inundation_ID')
    176250
    177251   ####################################################################################
    … …  
    193267##    if do_all!=True:
    194268##        print 'continuing'
    195     if gp.Exists(energy_less_final_dis):
     269    if (gp.Exists(energy_less_final_dis) and do_all == False):
    196270        print 'analysis already done, skipping to next one'
    197271    else:
    198272        # Process: Extract elevation data by Mask...
    199273        print 'check if elevation already extracted'
    200         if gp.Exists(elevation_extract)!=True:
     274        if (gp.Exists(elevation_extract)!=True or do_all == True):
    201275            print 'not already extracted'
    202276            print 'extracting by mask'
    … …  
    207281        # Process: Reclassify elevation data
    208282        print 'check if elevation already reclassified'
    209         if gp.Exists(energy_reclass)!=True:
     283        if (gp.Exists(energy_reclass)!=True or do_all == True):
    210284            print 'reclassify based on elevation height of:', energy_run_up_height
    211285            reclass_level = "-200 "+energy_run_up_height+" 1;"+energy_run_up_height+" 300 0"
    … …  
    216290        # Process: Raster to Polygon...
    217291        print 'check if already converted from raster to polyon'
    218         if gp.Exists(energy_poly)!=True:
     292        if (gp.Exists(energy_poly)!=True or do_all == True):
    219293            print 'raster to polyon'
    220294            gp.RasterToPolygon_conversion(energy_reclass, energy_poly, "NO_SIMPLIFY", "VALUE")
    … …  
    254328    inundation = raster_path + "\\"+ inundation_name[0]
    255329    inundation_area1 = "N:\\georisk_models\\inundation\data\\new_south_wales\\batemans_bay_tsunami_scenario_2009\\anuga\\polygons\\polygons.gdb\inundation_area1"
    256     inundation_extract = "inundation_"+run_up_name+"_extract"
    257     inundation_reclass = "inundation_"+run_up_name+"_reclass"
    258     inundation_poly = "inundation_"+run_up_name+"_poly"
    259     inundation_less = "inundation_"+run_up_name+"_less"
     330    inundation_extract = "inundation_extract"
     331    inundation_reclass = "inundation_reclass"
     332    inundation_poly = "inundation_poly"
     333    inundation_less = "inundation_less"
    260334    inundation_less_final = "inundation_less_final"
    261335    inundation_less_final_dis = "inundation_less_final_dis"
    262336
    263337    # Check if already done...
    264     if gp.Exists(inundation_less_final_dis):
     338    if (gp.Exists(inundation_less_final_dis) and do_all == False):
    265339        print 'analysis already done, skipping to next one'
    266340    else:
    267341        # Process: Extract inundation data by Mask...
    268342        print 'check if inundation already extracted'
    269         if gp.Exists(inundation_extract)!=True:
     343        if (gp.Exists(inundation_extract)!=True or do_all == True):
    270344            print 'not already extracted'
    271345            print 'extracting by mask'
    … …  
    276350        # Process: Reclassify inundation data
    277351        print 'check if inundation already reclassified'
    278         if gp.Exists(inundation_reclass)!=True:
     352        if (gp.Exists(inundation_reclass)!=True or do_all == True):
    279353            print 'reclassify based on inundation'
    280354            gp.Reclassify_sa(inundation_extract, "Value", "-200 0 1;0 300 0", inundation_reclass, "DATA")
    … …  
    284358        # Process: Raster to Polygon...
    285359        print 'check if already converted from raster to polyon'
    286         if gp.Exists(inundation_poly)!=True:
     360        if (gp.Exists(inundation_poly)!=True or do_all == True):
    287361            print 'raster to polyon'
    288362            gp.RasterToPolygon_conversion(inundation_reclass, inundation_poly, "NO_SIMPLIFY", "VALUE")
    … …  
    329403
    330404    # Search feature class for areas
     405    print 'get Madsen surf similarity areas'
     406    surf_sim_area_dict = {}
     407    #analytical_areas = []
     408    cur = gp.SearchCursor(elevation_surf_sim_less_final_dis)
     409    sRow = cur.Next()
     410    while sRow:
     411        if not(sRow.GetValue("inundation_ID")in surf_sim_area_dict):
     412            surf_sim_area_dict[sRow.GetValue("inundation_ID")] = sRow.GetValue("Shape_Area")
     413        else:
     414            surf_sim_area_dict[sRow.GetValue("inundation_ID")] = (surf_sim_area_dict[sRow.GetValue("inundation_ID")]+sRow.GetValue("Shape_Area"))
     415        #analytical_areas.append(sRow.GetValue("Shape_Area"))
     416        sRow = cur.Next()
     417    #print analytical_areas
     418    print surf_sim_area_dict
     419
     420    # Search feature class for areas
    331421    print 'get Energy areas'
    332422    energy_area_dict = {}
    … …  
    367457        if not key in analytical_area_dict:
    368458            #analytical_area_dict[key] = None
    369             area_comp_madsen[key].append(None)
    370             area_percent_madsen[key].append(None)
     459            area_comp_madsen[key].append(-999999)
     460            area_percent_madsen[key].append(-999999)
    371461        else:
    372462            diff_madsen = anuga_area_dict[key]-analytical_area_dict[key]
    … …  
    374464            area_comp_madsen[key].append(diff_madsen)
    375465            area_percent_madsen[key].append(percent_diff_madsen)
     466
     467        if not key in surf_sim_area_dict:
     468            #analytical_area_dict[key] = None
     469            area_comp_surf_sim[key].append(-999999)
     470            area_percent_surf_sim[key].append(-999999)
     471        else:
     472            diff_surf_sim = anuga_area_dict[key]-surf_sim_area_dict[key]
     473            percent_diff_surf_sim = 100*diff_surf_sim/anuga_area_dict[key]
     474            area_comp_surf_sim[key].append(diff_surf_sim)
     475            area_percent_surf_sim[key].append(percent_diff_surf_sim)
    376476           
    377477        if not key in energy_area_dict:
    378             area_comp_energy[key].append(None)
    379             area_percent_energy[key].append(None)
     478            area_comp_energy[key].append(-999999)
     479            area_percent_energy[key].append(-999999)
    380480        else:
    381481            diff_energy = anuga_area_dict[key]-energy_area_dict[key]
    … …  
    392492csv_anuga_sum = csv.writer(open(out_file_anuga_sum,'w'))
    393493csv_madsen = csv.writer(open(out_file_madsen,'w'))
     494csv_surf_sim = csv.writer(open(out_file_surf_sim,'w'))
    394495csv_energy = csv.writer(open(out_file_energy,'w'))
    395496
    396 csv_anuga.writerow([1,2,3,4,5,6,7,8,9])
    397 csv_madsen.writerow([1,2,3,4,5,6,7,8,9])
    398 csv_energy.writerow([1,2,3,4,5,6,7,8,9])
     497csv_anuga.writerow([1,2,3,4,5,6,7,8,9,10])
     498csv_madsen.writerow([1,2,3,4,5,6,7,8,9,10])
     499csv_energy.writerow([1,2,3,4,5,6,7,8,9,10])
     500csv_surf_sim.writerow([1,2,3,4,5,6,7,8,9,10])
    399501
    400502for key, value in anuga_sum_areas.iteritems():
    … …  
    408510            area_total_anuga[key][i]
    409511        except NameError:
    410             value_list_anuga.append(None)
     512            value_list_anuga.append(-9999)
    411513        else:
    412514            value_list_anuga.append(area_total_anuga[key][i]) 
    … …  
    419521            area_percent_madsen[key][i]
    420522        except NameError:
    421             value_list_madsen.append(None)
     523            value_list_madsen.append(-9999)
    422524        else:
    423525            value_list_madsen.append(area_percent_madsen[key][i])
    424526    csv_madsen.writerow(value_list_madsen)
     527
     528for i in range(len(area_percent_surf_sim[1])):
     529    value_list_surf_sim =[]
     530    for key in area_percent_surf_sim:
     531        try:
     532            area_percent_surf_sim[key][i]
     533        except NameError:
     534            value_list_surf_sim.append(-9999)
     535        else:
     536            value_list_surf_sim.append(area_percent_surf_sim[key][i])
     537    csv_surf_sim.writerow(value_list_surf_sim)
    425538   
    426539for i in range(len(area_percent_energy[1])):
    … …  
    430543            area_percent_energy[key][i]
    431544        except NameError:
    432             value_list_energy.append(None)
     545            value_list_energy.append(-9999)
    433546        else:
    434547            value_list_energy.append(area_percent_energy[key][i]) 

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

    r7577 r7592  
    99boundary_path = "/nas/gemd/georisk_models/inundation/data/new_south_wales/batemans_bay_tsunami_scenario_2009/anuga/boundaries/"
    1010in_file_name = 'area_comparisons_madsen.csv'
    11 #in_file_name = 'area_comparisons_energy.csv'
     11##in_file_name = 'area_comparisons_energy.csv'
     12##in_file_name = 'area_comparisons_surf_sim.csv'
    1213figure_name = 'area_comparisons_madsen.png'
     14##figure_name = 'area_comparisons_energy.png'
     15##figure_name = 'area_comparisons_surf_sim.png'
    1316
    1417in_file = join(boundary_path, in_file_name)
    … …  
    3841##        key_list.append(key)
    3942pylab.scatter(area_id_list,area_percent_list)
    40 pylab.axis([0,10,-2000,2000])
     43pylab.axis([0,11,-500,500])
    4144pylab.savefig(figure_path)

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

    r7577 r7592  
    3636              'Event3_MSL':58284,
    3737              'Puysegur_200yr':58129,
    38               'Puysegur_500yr':58115,
     38              #'Puysegur_500yr':58115,
    3939              'Puysegur_1000yr':58226,
    4040              'Puysegur_5000yr':58286,

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

    r7577 r7592  
    88from time import localtime, strftime, gmtime
    99from os.path import join, exists
     10import sys
    1011
    1112
    … …  
    2627# Model specific parameters.
    2728# One or all can be changed each time the run_model script is executed
    28 tide = 1.0                # difference between MSL and HAT (1.0)
     29tide = 0.0                # difference between MSL and HAT (1.0)
    2930
    3031# the event number or the mux file name
    … …  
    3435##event_number = 58284     #1 in 2000 yr Puysegur (Event 3)
    3536##event_number = 58286     #1 in 5000 yr Puysegur
    36 event_number = 58346     #1 in 10000 yr Puysegur (Event 1)
     37##event_number = 58346     #1 in 10000 yr Puysegur (Event 1)
    3738
    3839##event_number = 51077     #1 in 200 yr New Hebrides
    … …  
    4748######event_number = 58272      #1 in ~15000 yr Puysegur
    4849######event_number = 51445      #1 in ~15000 yr New Hebrides
     50
     51if len(sys.argv) > 1:
     52    event_number = int(sys.argv[1])
     53else:   
     54    event_number = 58346
     55##event_number_list = [7871,31583,31602,31714,31913] #1.5-2.0m
     56##event_number_list = [31997,51268,51375,51436,51452] #1.5-2.0m
     57##event_number_list =  [58176,58286,58318,58337,58350] #1.5-2.0m
     58
     59##event_number_list = [31708,31831,31897,31961,31973]
     60##event_number_list = [31981,51270,51328,51415,51460]
     61event_number_list = [51470,58274,58331,58348,58357,58360]
     62
    4963
    5064alpha = 0.1             # smoothing parameter for mesh

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

    r7586 r7592  
    1414
    1515
    16 def wave_energy(filepath, filebasename):
     16def wave_energy(filepath, filebasename, slope):
    1717   
    1818    filepathlist = []
    … …  
    152152        # call run-up module
    153153        #####################################################################
    154         Run_up = run_up.Madsen(100,0.017,max_stage,max_time)
     154        Run_up = run_up.Madsen(100,numpy.tan(slope),max_stage,max_time)
    155155        print 'Madsen runup', Run_up
    156156        run_up_list.append(Run_up)
    157157
    158         Run_up_surf_sim = run_up.surf_sim(100,0.03,max_stage,max_time)
     158        Run_up_surf_sim = run_up.surf_sim(100,numpy.tan(slope),max_stage,max_time)
    159159        print 'Madsen surf similarity runup', Run_up_surf_sim
    160160        run_up_list_surf_sim.append(Run_up_surf_sim)
    161161       
    162         energy_run_up  = run_up.energy_conserv(max_energy,1.0,0.8,9.81)
     162        energy_run_up  = run_up.energy_conserv(max_energy,slope+0.2,slope,9.81)
    163163        energy_run_up_list.append(energy_run_up)
    164164        print 'energy run up', energy_run_up
    … …  
    182182    outFilename= filepath +'/wave_energy_summary.csv'
    183183    outFile = file(outFilename,'w')
    184     outFile.write('filename, max crest-integrated energy (J.s/m^2), max instantatneous energy J/m^2, Run up (Madsen), Energy run up , max stage (m)\n')
     184    outFile.write('filename, max crest-integrated energy (J.s/m^2), max instantatneous energy J/m^2, Run up (Madsen), Energy run up , max stage (m), Surf similarity run up (m)\n')
    185185    for filename in filepathlist:
    186186        outFile.write(filename+',')
    … …  
    229229    index = 2872
    230230    filebasename = 'sts_gauge_' + str(index) +'.csv'
    231     wave_energy(filepath, filebasename)
     231    slope = 1 #degrees
     232    wave_energy(filepath, filebasename, slope)
Note: See TracChangeset for help on using the changeset viewer.