Changeset 7570


Ignore:
Timestamp:
Nov 25, 2009, 4:28:13 PM (9 years ago)
Author:
jgriffin
Message:
 
File:
1 edited

Legend:

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

    r7568 r7570  
     1"""# ---------------------------------------------------------------------------
     2 elevation_model.py
     3 Created on: Tue Nov 17 2009 02:11:08
     4 Created by: Jonathan Griffin
     5 This script compares calculated inundate areas for run-up solutions from analystical models
     6 and ANUGA modelling.
     7
     8
    19# ---------------------------------------------------------------------------
    2 # elevation_model.py
    3 # Created on: Tue Nov 17 2009 02:11:08
    4 #   (generated by ArcGIS/ModelBuilder)
    5 # ---------------------------------------------------------------------------
    6 
     10"""
    711# Import system modules
    812import sys, string, os, arcgisscripting
    913from os.path import join
     14from collections import defaultdict
    1015# Create the Geoprocessor object
    1116gp = arcgisscripting.create(9.3)
     
    4348##              'New_Hebrides_2000yr':51292,
    4449##              'New_Hebrides_5000yr':51424,}
     50
     51# Dictionary for storing area comparison
     52area_comp = defaultdict(list)
    4553
    4654boundaries = file(boundary_file,'r')
     
    7280    ##run_up_madsen = 6.0
    7381    run_up_height = str(run_up_madsen)
    74     run_up_name = run_up_height.split('.')[0]
     82    run_up_name = run_up_height.replace('.','_')[0]
    7583
    7684
     
    178186    # Process: Select...
    179187    print 'select by attribute'
    180     gp.Select_analysis(inundation_poly, inundation_less, "\"GRIDCODE\" = 1")
     188    gp.Select_analysis(inundation_poly, inundation_less, "\"GRIDCODE\" = 0")
    181189
    182190    # Process: Create layer...
     
    197205    # Search feature class for areas
    198206    print 'get areas'
    199     area_dict = {}
    200     areas = []
     207    analytical_area_dict = {}
     208    #analytical_areas = []
    201209    cur = gp.SearchCursor(elevation_less_final)
    202210    sRow = cur.Next()
    203211    while sRow:
    204         area_dict[sRow.GetValue("inundation_ID")] = sRow.GetValue("Shape_Area")
    205         areas.append(sRow.GetValue("Shape_Area"))
     212        if not(sRow.GetValue("inundation_ID")in analytical_area_dict):
     213            analytical_area_dict[sRow.GetValue("inundation_ID")] = sRow.GetValue("Shape_Area")
     214        else:
     215            analytical_area_dict[sRow.GetValue("inundation_ID")] = (analytical_area_dict[sRow.GetValue("inundation_ID")]+sRow.GetValue("Shape_Area"))
     216        #analytical_areas.append(sRow.GetValue("Shape_Area"))
    206217        sRow = cur.Next()
    207     print areas
    208     print area_dict
    209 
     218    #print analytical_areas
     219    print analytical_area_dict
     220
     221    # Search feature class for areas
     222    print 'get areas'
     223    anuga_area_dict = {}
     224    #anuga_areas = []
     225    cur = gp.SearchCursor(inundation_less_final)
     226    sRow = cur.Next()
     227    while sRow:
     228        if not(sRow.GetValue("inundation_ID") in anuga_area_dict):
     229            anuga_area_dict[sRow.GetValue("inundation_ID")]=sRow.GetValue("Shape_Area")
     230        else:
     231            anuga_area_dict[sRow.GetValue("inundation_ID")]= anuga_area_dict[sRow.GetValue("inundation_ID")]+sRow.GetValue("Shape_Area") 
     232        #anuga_areas.append(sRow.GetValue("Shape_Area"))
     233        sRow = cur.Next()
     234    #print anuga_areas
     235    print anuga_area_dict
     236
     237    for key in anuga_area_dict:
     238        diff = anuga_area_dict[key]-analytical_area_dict[key]
     239        area_comp[key].append(diff)
     240    print 'area_comp', area_comp
     241   
    210242boundaries.close()
Note: See TracChangeset for help on using the changeset viewer.