Changeset 3721


Ignore:
Timestamp:
Oct 10, 2006, 11:36:20 AM (17 years ago)
Author:
sexton
Message:

hobart testing and report making

Location:
anuga_work/production/hobart_2006
Files:
16 added
4 edited

Legend:

Unmodified
Added
Removed
  • anuga_work/production/hobart_2006/export_results.py

    r3680 r3721  
    66from os import sep
    77
    8 # OLE - this is the smallest sww file with the complex bounding polygon
    9 time_dir = '20060928_073318' #alpha = 0.5
     8time_dir = '20061006_062319' # Mw 8-5
    109
    11 #time_dir = '20060928_064732' #alpha = 0.01
    12 #time_dir = '20060928_220905' #alpha = 1.0
    13 #time_dir = '20060926_064750' #this worked for cellsize 30 and 100
    14 #time_dir = '20060929_033100' #alpha = 0.01 poly anticlockwise - DIDN'T WORK
    15 #time_dir = '20060929_075954' #alpha = 0.1 grid with MOST
    16 #time_dir = '20060929_075009' #alpha = 0.1 points with MOST
    1710directory = project.outputdir
    1811name = directory + time_dir +sep + 'source'
    1912
    2013print 'output dir:', name
    21 #print 'Which variable do you want to export?'
    22 #which_var = int(raw_input('Stage = 0, Absolute Momentum = 1, Depth = 2, Speed = 3  '  ))
    2314which_var = 4
    2415if which_var == 0:  # Stage
     
    3122
    3223if which_var == 2:  # Depth
    33     outname = name + '_depth'
     24    outname = name + '_depth_bruny'
    3425    quantityname = 'stage-elevation'  #Depth
    3526
     
    3930
    4031if which_var == 4:  # Elevation
    41     outname = name + '_elevation_2'
     32    outname = name + '_elevation'
    4233    quantityname = 'elevation'  #Elevation
    4334
    4435print 'start sww2dem'
     36# for hobart
     37##sww2dem(name, basename_out = outname,
     38##            quantity = quantityname,
     39##            cellsize = 25,      # would prefer this at 25
     40##            # define region for viz purposes
     41##            easting_min = project.eastingmin25,
     42##            easting_max = project.eastingmax25,
     43##            northing_min = project.northingmin25,
     44##            northing_max = project.northingmax25,       
     45##            reduction = max, #this is because we want max quantityname
     46##            verbose = True,
     47##            format = 'asc')
     48
     49# for bruny
    4550sww2dem(name, basename_out = outname,
    4651            quantity = quantityname,
    47             cellsize = 500,      # would prefer this at 25
     52            cellsize = 25,      # would prefer this at 25
    4853            # define region for viz purposes
    49             #easting_min = project.e_min_area,
    50             #easting_max = project.e_max_area,
    51             #northing_min = project.n_min_area,
    52             #northing_max = project.n_max_area,       
     54            easting_min = project.eastingmin25_2,
     55            easting_max = project.eastingmax25_2,
     56            northing_min = project.northingmin25_2,
     57            northing_max = project.northingmax25_2,       
    5358            reduction = max, #this is because we want max quantityname
    5459            verbose = True,
    5560            format = 'asc')
    56 
  • anuga_work/production/hobart_2006/get_timeseries.py

    r3687 r3721  
    2323    mkdir (reportdir)
    2424   
    25 # data
    26 production_dirs = {'20061003_062209': 'test'}
    27                    #'20061003_062209': 'Mw 8.7'}#,
    28                    #'': 'Mw 8.5'}
     25# 2500 interior resolution for both events
     26# use both to compare the events
     27production_dirs = {'20061006_062319': 'Mw 8.5'}
     28                   #'20061008_234702': 'Mw 8.7'}
    2929
    3030# Generate figures
  • anuga_work/production/hobart_2006/project.py

    r3701 r3721  
    77import sys
    88from time import localtime, strftime, gmtime
    9 from anuga.utilities.polygon import read_polygon, plot_polygons
     9from anuga.utilities.polygon import read_polygon, plot_polygons, polygon_area
    1010from anuga.coordinate_transforms.redfearn import degminsec2decimal_degrees, convert_points_from_latlon_to_utm
    1111
     
    2626scenario_dir_name = 'hobart_tsunami_scenario_2006'
    2727
     28# data provided by Tas SES and checked by NM&I
     29#onshore_name = 'hob3_topo' # original
     30onshore_name_25 = 'hob5_topo_25m' # 25m grid and clipped to 100m elevation or 3000m from coast
     31onshore_name = 'hob5_topo' # 12.5m grid and clipped to 100m elevation or 3000m from coast
     32#onshore_name_25 = 'hob6_topo_25m' # 25m grid NOT clipped
     33offshore_name_tas1 = 'derwent_2m'
     34offshore_name_tas2 = 'derwent_5m'
     35offshore_name_tas3 = 'south_east_tas' #actually this is AHO
     36offshore_name_tas4 = 'hobart_1m'
    2837
    29 ##
    30 ### data provided by Tas SES and checked by NM&I
    31 ###onshore_name = 'hob3_topo' # original
    32 ###onshore_name_25 = 'hob5_topo_25m' # 25m grid and clipped to 100m elevation or 3000m from coast
    33 ##onshore_name = 'hob5_topo' # 12.5m grid and clipped to 100m elevation or 3000m from coast
    34 ##onshore_name_25 = 'hob6_topo_25m' # 25m grid NOT clipped
    35 ##offshore_name_tas1 = 'derwent_2m'
    36 ##offshore_name_tas2 = 'derwent_5m'
    37 ##offshore_name_tas3 = 'south_east_tas' #actually this is AHO
    38 ##offshore_name_tas4 = 'hobart_1m'
    39 ##
    40 ### AHO data and checked by NM&I
    41 ##offshore_name1 = 'xy100003760'
    42 ##offshore_name2 = 'xy100003761'
    43 ##offshore_name3 = 'xy100003762'
    44 ##offshore_name4 = 'xy100003907'
    45 ##offshore_name5 = 'xy100003908'
    46 ##offshore_name6 = 'xy100003909'
    47 ##offshore_name7 = 'xy100003910'
    48 ##offshore_name8 = 'xy100003932'
    49 ##offshore_name9 = 'xy100003933'
    50 ##offshore_name10 = 'xy100003934'
    51 ##offshore_name11 = 'xy100003935'
    52 ##offshore_name12 = 'xy100003936'
    53 ##offshore_name13 = 'xy100003964'
    54 ##offshore_name14 = 'xy100014250'
    55 ##offshore_name15 = 'xy100014253'
    56 ##offshore_name16 = 'xy100016142'
    57 ##
    58 ### developed by NM&I
    59 ##coast_name = 'coastline_points'
     38# AHO data and checked by NM&I
     39offshore_name1 = 'xy100003760'
     40offshore_name2 = 'xy100003761'
     41offshore_name3 = 'xy100003762'
     42offshore_name4 = 'xy100003907'
     43offshore_name5 = 'xy100003908'
     44offshore_name6 = 'xy100003909'
     45offshore_name7 = 'xy100003910'
     46offshore_name8 = 'xy100003932'
     47offshore_name9 = 'xy100003933'
     48offshore_name10 = 'xy100003934'
     49offshore_name11 = 'xy100003935'
     50offshore_name12 = 'xy100003936'
     51offshore_name13 = 'xy100003964'
     52offshore_name14 = 'xy100014250'
     53offshore_name15 = 'xy100014253'
     54offshore_name16 = 'xy100016142'
    6055
    61 boundary_basename = 'puysegur' # Mw 8.7
    62 #boundary_basename = 'puysegur_clip' # Mw 8.5
     56# developed by NM&I
     57coast_name = 'coastline_points'
    6358
     59#boundary_basename = 'puysegur' # Mw 8.7
     60boundary_basename = 'puysegur_clip' # Mw 8.5
    6461
    6562#swollen/ all data output
     
    8178
    8279gauge_filename = gaugedir + 'hobart_gauges_final.csv'
    83 #buildings_filename = gaugedir + 'hobart_res.csv'
    84 #buildings_filename_damage_out = 'hobart_res_modified.csv'
    85 #gaugetimeseries = gaugedir + 'hobart'
    86 
    87 # boundary source data
    88 #MOST_dir = 'f:'+sep+'3'+sep+'ehn'+sep+'users'+sep+'davidb'+sep+'tsunami'+sep+'WA_project'+sep+'SU-AU_90'+sep+'most_2'+sep+'detailed'+sep
    8980
    9081codedir = getcwd()+sep                           
     
    9384
    9485# Necessary if using point datasets, rather than grid
    95 ##onshore_dem_name = datadir + onshore_name
    96 ##onshore_dem_name_25 = datadir + onshore_name_25
    97 ##all_onshore_dem_name = datadir + 'combined_onshore'
    98 ##offshore_dem_name_local1 = datadir + offshore_name_tas1
    99 ##offshore_dem_name_local2 = datadir + offshore_name_tas2
    100 ##offshore_dem_name_local3 = datadir + offshore_name_tas3
    101 ##offshore_dem_name_local4 = datadir + offshore_name_tas4
    102 ##offshore_dem_name_aho1 = datadir + offshore_name1
    103 ##offshore_dem_name_aho2 = datadir + offshore_name2
    104 ##offshore_dem_name_aho3 = datadir + offshore_name3
    105 ##offshore_dem_name_aho4 = datadir + offshore_name4
    106 ##offshore_dem_name_aho5 = datadir + offshore_name5
    107 ##offshore_dem_name_aho6 = datadir + offshore_name6
    108 ##offshore_dem_name_aho7 = datadir + offshore_name7
    109 ##offshore_dem_name_aho8 = datadir + offshore_name8
    110 ##offshore_dem_name_aho9 = datadir + offshore_name9
    111 ##offshore_dem_name_aho10 = datadir + offshore_name10
    112 ##offshore_dem_name_aho11 = datadir + offshore_name11
    113 ##offshore_dem_name_aho12 = datadir + offshore_name12
    114 ##offshore_dem_name_aho13 = datadir + offshore_name13
    115 ##offshore_dem_name_aho14 = datadir + offshore_name14
    116 ##offshore_dem_name_aho15 = datadir + offshore_name15
    117 ##offshore_dem_name_aho16 = datadir + offshore_name16
    118 ##coast_dem_name = datadir + coast_name
    119 
     86onshore_dem_name = datadir + onshore_name
     87onshore_dem_name_25 = datadir + onshore_name_25
     88all_onshore_dem_name = datadir + 'combined_onshore'
     89offshore_dem_name_local1 = datadir + offshore_name_tas1
     90offshore_dem_name_local2 = datadir + offshore_name_tas2
     91offshore_dem_name_local3 = datadir + offshore_name_tas3
     92offshore_dem_name_local4 = datadir + offshore_name_tas4
     93offshore_dem_name_aho1 = datadir + offshore_name1
     94offshore_dem_name_aho2 = datadir + offshore_name2
     95offshore_dem_name_aho3 = datadir + offshore_name3
     96offshore_dem_name_aho4 = datadir + offshore_name4
     97offshore_dem_name_aho5 = datadir + offshore_name5
     98offshore_dem_name_aho6 = datadir + offshore_name6
     99offshore_dem_name_aho7 = datadir + offshore_name7
     100offshore_dem_name_aho8 = datadir + offshore_name8
     101offshore_dem_name_aho9 = datadir + offshore_name9
     102offshore_dem_name_aho10 = datadir + offshore_name10
     103offshore_dem_name_aho11 = datadir + offshore_name11
     104offshore_dem_name_aho12 = datadir + offshore_name12
     105offshore_dem_name_aho13 = datadir + offshore_name13
     106offshore_dem_name_aho14 = datadir + offshore_name14
     107offshore_dem_name_aho15 = datadir + offshore_name15
     108offshore_dem_name_aho16 = datadir + offshore_name16
     109coast_dem_name = datadir + coast_name
    120110
    121111# addition once total grid delivered
     
    128118combined_dem_name   = datadir + 'hobart_combined_elevation'
    129119combined_dem_name_2 = datadir + 'hobart_combined_elevation_2'
     120combined_dem_name_3 = datadir + 'hobart_combined_elevation_jane'
     121combined_dem_name_4 = datadir + 'hobart_combined_elevation_test'
    130122
    131123#outputname = outputtimedir + basename  #Used by post processing
     
    150142refzone = zone
    151143
    152 # Second run - bottom bright, topr, top, left
    153 #polyAll = [[520000, 5170000],[580000, 5170000],[580000, 5200000],[590000,5240000],[520000,5260000]]
    154 
    155 
    156144# Mark run - morning Fri 29 Sep; surrounds -100m and 20mish elevation
    157 polyAll = read_polygon(polygondir+'new_extent.csv')
     145polyAll = read_polygon(polygondir+'new_extent_2.csv')
    158146plot_polygons([polyAll, polyAll2],'boundingpoly',verbose=False)
    159 
     147print 'Area of bounding polygon', polygon_area(polyAll)
    160148
    161149###################################################################
     
    163151###################################################################
    164152
     153# clipping 12.5m onshore data set
     154eastingmin = 520000
     155eastingmax = 536000
     156northingmin = 5245000
     157northingmax = 5260000
    165158
    166159# clipping 25m data set - Hobart
     
    176169northingmax25_2 = 5212052.309
    177170
    178 
    179 
    180171###############################
    181172# Interior region definitions
    182173###############################
    183 
    184 ###Interior region - Hobart city area + Glenorchy, Kingston
    185 ##i0 = [517000, 5267000]
    186 ##i1 = [517000, 5255000]
    187 ##i2 = [520000, 5250000]
    188 ##i3 = [522000, 5239000]
    189 ##i4 = [524000, 5238000]
    190 ##i5 = [526000, 5236000]
    191 ##i6 = [530000, 5244000]
    192 ##i7 = [530000, 5250000]
    193 ##i8 = [534000, 5254000]
    194 ##i9 = [520000, 5270000]
    195 ##
    196 ##poly_hobart = [i0, i1, i2, i3, i4, i5, i6, i7, i8, i9]
    197 ##
    198 ### Tasman Peninsula
    199 ##l0 = [550000, 5247000]
    200 ##l1 = [550000, 5211000]
    201 ##l2 = [583000, 5211000]
    202 ##l3 = [583000, 5247000]
    203 ##
    204 ##poly_tasman_peninsula = [l0, l1, l2, l3]
    205 ##
    206 ### Bruny Island
    207 ##poly_bruny = read_polygon(polygondir+'bruny.csv')
    208 
    209174
    210175# Hobart digitized polygons
     
    212177poly_hobart2 = read_polygon(polygondir+'Hob_poly2.csv')
    213178poly_hobart3 = read_polygon(polygondir+'Hob_poly3.csv')
     179poly_hobart4 = read_polygon(polygondir+'Hob_poly4.csv')
    214180
    215 plot_polygons([polyAll, poly_hobart1,poly_hobart2,poly_hobart3],'boundingpoly2',verbose=False)
     181plot_polygons([polyAll, poly_hobart1,poly_hobart2,poly_hobart3,poly_hobart4],'boundingpoly2',verbose=False)
  • anuga_work/production/hobart_2006/run_hobart.py

    r3701 r3721  
    7373# creates pts file for combined 50m DEM
    7474dem2pts(onshore_offshore_dem_name, use_cache=True, verbose=True)
    75 
     75"""
    7676# 25m data (clipping the around the Hobart area)
    7777convert_dem_from_ascii2netcdf(onshore_offshore_dem_name_25, use_cache=True, verbose=True)
     
    113113#G.export_points_file(project.combined_dem_name + '.pts')
    114114
    115 
     115"""
    116116#----------------------------------------------------------------------------
    117117# Create the triangular mesh based on overall clipping polygon with a tagged
     
    146146                              'e12': [12], 'e13': [13], 'e14': [14],
    147147                              'e15': [15]},
    148            'maximum_triangle_area': 750000,
    149            'filename': meshname,
    150            'interior_regions': interior_regions},
     148           'maximum_triangle_area': 2500000,
     149           'filename': meshname},
     150           #'interior_regions': interior_regions},
    151151          verbose = True, evaluate=False)
    152152
     
    176176domain.set_quantity('friction', 0.0)
    177177
    178 domain.set_quantity('elevation',
    179                     filename = project.combined_dem_name_2 + '.pts',
    180                     use_cache = True,
    181                     verbose = True,
    182                     alpha = 0.1
    183                     )
     178#domain.set_quantity('elevation',
     179#                    filename = project.combined_dem_name_2 + '.pts',
     180#                    filename = onshore_offshore_dem_name + '.pts',
     181#                    use_cache = True,
     182#                    verbose = True,
     183#                    alpha = 0.1
     184#                    )
    184185
    185186
     
    188189#-------------------------------------------------------------------------------
    189190print 'start ferret2sww'
     191print '', project.boundary_basename
    190192from anuga.shallow_water.data_manager import ferret2sww
    191193
     
    195197east  = project.east
    196198
    197 ###note only need to do when an SWW file for the MOST boundary doesn't exist
    198 ##cache(ferret2sww,
    199 ##      (source_dir + project.boundary_basename,
    200 ##       source_dir + project.boundary_basename),
    201 ##      {'verbose': True,
    202 ##       'minlat': south,
    203 ##       'maxlat': north,
    204 ##       'minlon': west,
    205 ##       'maxlon': east,
    206 ###       'origin': project.mesh_origin,
    207 ##       'origin': domain.geo_reference.get_origin(),
    208 ##       'mean_stage': tide,
    209 ##       'zscale': 1,                 #Enhance tsunami
    210 ##       'fail_on_NaN': False,
    211 ##       'inverted_bathymetry': True},
    212 ##      #evaluate = True,
    213 ##       verbose = True,
    214 ##       dependencies = source_dir + project.boundary_basename + '.sww')
     199#note only need to do when an SWW file for the MOST boundary doesn't exist
     200cache(ferret2sww,
     201      (source_dir + project.boundary_basename,
     202       source_dir + project.boundary_basename),
     203      {'verbose': True,
     204       'minlat': south,
     205       'maxlat': north,
     206       'minlon': west,
     207       'maxlon': east,
     208#       'origin': project.mesh_origin,
     209       'origin': domain.geo_reference.get_origin(),
     210       'mean_stage': tide,
     211       'zscale': 1,                 #Enhance tsunami
     212       'fail_on_NaN': False,
     213       'inverted_bathymetry': True},
     214      #evaluate = True,
     215       verbose = True,
     216       dependencies = source_dir + project.boundary_basename + '.sww')
    215217
    216218
Note: See TracChangeset for help on using the changeset viewer.