Changeset 3671


Ignore:
Timestamp:
Sep 27, 2006, 5:55:02 PM (18 years ago)
Author:
sexton
Message:

Hobart updates; using gridded data

Location:
anuga_work/production/hobart_2006
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • anuga_work/production/hobart_2006/export_results.py

    r3668 r3671  
    33
    44from anuga.shallow_water.data_manager import sww2dem
    5 #from anuga.pyvolution.ermapper_grids import read_ermapper_grid
    65from anuga.abstract_2d_finite_volumes.util import Screen_Catcher
    76from os import sep
    87
    9 time_dir = '20060925_115139' #MSL
     8time_dir = '20060927_054931' #MSL
    109directory = project.outputdir
    1110name = directory + time_dir +sep + 'source'
    … …  
    1514#which_var = int(raw_input('Stage = 0, Absolute Momentum = 1, Depth = 2, Speed = 3  '  ))
    1615which_var = 4
    17 #sys.stderr.write(sys.stdout.data)
    1816if which_var == 0:  # Stage
    1917    outname = name + '_stage'
    … …  
    3028if which_var == 3:  # Speed
    3129    outname = name + '_speed'
    32     #quantityname = '((xmomentum/(stage-elevation))**2 + (ymomentum/(stage-elevation))**2)**0.5'  #Speed
    33     #quantityname = 'xmomentum/(stage-elevation)'  #Speed
    3430    quantityname = '(xmomentum**2 + ymomentum**2)**0.5/(stage-elevation+1.e-30)'  #Speed
    3531
    3632if which_var == 4:  # Elevation
    37     outname = name + '_elevation_node3'
     33    outname = name + '_elevation'
    3834    quantityname = 'elevation'  #Elevation
    3935
    4036print 'start sww2dem'
    41 #sys.stderr.write(sys.stdout.data)
    4237sww2dem(name, basename_out = outname,
    4338            quantity = quantityname,
    44             cellsize = 30,      # would prefer this at 25
     39            cellsize = 25,      # would prefer this at 25
    4540            # define region for viz purposes
    4641            #easting_min = project.e_min_area,
    … …  
    5247            format = 'asc')
    5348
    54 #sys.stderr.write(sys.stdout.data)
    55 
    56 #Check
    57 
    58 #data = read_ermapper_grid(name)
    59 #print 'Values from %s are in [%f, %f]' %(name, min(data.flat), max(data.flat))

  • anuga_work/production/hobart_2006/project.py

    r3669 r3671  
    77import sys
    88from time import localtime, strftime, gmtime
    9                
     9from anuga.utilities.polygon import read_polygon, plot_polygons
     10from anuga.coordinate_transforms.redfearn import degminsec2decimal_degrees, convert_points_from_latlon_to_utm
     11
     12if sys.platform == 'win32':
     13    home = getenv('INUNDATIONHOME')
     14    user = getenv('USERPROFILE')
     15
     16else:   
     17    home = getenv('INUNDATIONHOME', sep+'d'+sep+'xrd'+sep+'gem'+sep+'2'+sep+'ramp'+sep+'risk_assessment_methods_project'+sep+'inundation')     
     18    user = getenv('LOGNAME')
     19    print 'USER:', user
     20
     21# INUNDATIONHOME is the inundation directory, not the data directory.
     22home += sep +'data'
     23
    1024#Making assumptions about the location of scenario data
    1125state = 'tasmania'
    … …  
    5064codename = 'project.py'
    5165
    52 if sys.platform == 'win32':
    53     home = getenv('INUNDATIONHOME')
    54     user = getenv('USERPROFILE')
    55 
    56 else:   
    57     home = getenv('INUNDATIONHOME', sep+'d'+sep+'xrd'+sep+'gem'+sep+'2'+sep+'ramp'+sep+'risk_assessment_methods_project'+sep+'inundation')     
    58     user = getenv('LOGNAME')
    59     print 'USER:', user
    60 
    61 # INUNDATIONHOME is the inundation directory, not the data directory.
    62 home += sep +'data'
    63 
    6466#Derive subdirectories and filenames
    6567#time = strftime('%Y%m%d_%H%M%S',localtime()) #gets time for new dir
    6668local_time = strftime('%Y%m%d_%H%M%S',gmtime()) #gets time for new dir
    67 print 'home', home
    6869meshdir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'meshes'+sep
    6970datadir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'topographies'+sep
    … …  
    7172polygondir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'polygons'+sep
    7273boundarydir = home+sep+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'boundaries'+sep
    73 #output dir without time
    7474outputdir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'outputs'+sep
    7575outputtimedir = outputdir + local_time + sep
    … …  
    7979buildings_filename = gaugedir + 'hobart_res.csv'
    8080buildings_filename_damage_out = 'hobart_res_modified.csv'
    81 
    8281gaugetimeseries = gaugedir + 'hobart'
    8382
    … …  
    8584#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
    8685
    87 codedir = getcwd()+sep
    88                                
     86codedir = getcwd()+sep                           
    8987codedirname = codedir + 'project.py'
    90 
    9188meshname = outputtimedir + 'mesh_' + basename
    9289
    … …  
    115112offshore_dem_name_aho16 = datadir + offshore_name16
    116113coast_dem_name = datadir + coast_name
     114# addition once total grid delivered
     115#onshore_offshore_dem_name = datadir + '25m_se_tas' #25m grid
     116onshore_offshore_dem_name = datadir + '50m_se_tas' #50m grid
     117
    117118combined_dem_name = datadir + 'hobart_combined_elevation'
     119combined_dem_name2 = datadir + 'hobart_combined_elevation2' # for alpha checking
    118120
    119121outputname = outputtimedir + basename  #Used by post processing
     122
     123###############################
     124# Domain definitions
     125###############################
     126
     127# bounding box
     128south = degminsec2decimal_degrees(-43,45,0)
     129north = degminsec2decimal_degrees(-42,30,0)
     130west = degminsec2decimal_degrees(146,45,0)
     131east = degminsec2decimal_degrees(148,0,0)#degminsec2decimal_degrees(148,15,0)
     132
     133#Main Domain of Hobart: first run JS 18/9/06
     134d0 = [south, west]
     135d1 = [south, east]
     136d2 = [north, east]
     137d3 = [north, west]
     138polyAll, zone = convert_points_from_latlon_to_utm([d0, d1, d2, d3])
     139refzone = zone
     140
     141# Second run - bottom bright, topr, top, left
     142polyAll = [[520000, 5170000],[580000, 5170000],[580000, 5200000],[590000,5240000],[520000,5260000]]
     143
     144# Third run - afternoon Wed 27 Sep; surrounds -100m and 20mish elevation
     145polyAll = read_polygon(polygondir+'new_extent.csv')
     146plot_polygons([polyAll],'boundingpoly',verbose=False)
     147
     148###################################################################
     149# Clipping regions for export to asc and regions for clipping data
     150###################################################################
     151
     152# region to export for inundation map
     153e_min_area = 521000#490000
     154e_max_area = 522000#580000
     155n_min_area = 5190000#5160000
     156n_max_area = 5200000#5275000
    120157
    121158# clipping 12.5m onshore data set
    … …  
    130167northingmax25 = 5260000
    131168
    132 # region to export for inundation map
    133 e_min_area = 521000#490000
    134 e_max_area = 522000#580000
    135 n_min_area = 5190000#5160000
    136 n_max_area = 5200000#5275000
    137 
    138 from anuga.coordinate_transforms.redfearn import degminsec2decimal_degrees, convert_points_from_latlon_to_utm
    139                        
    140 # bounding box
    141 south = degminsec2decimal_degrees(-43,45,0)
    142 north = degminsec2decimal_degrees(-42,30,0)
    143 west = degminsec2decimal_degrees(146,45,0)
    144 east = degminsec2decimal_degrees(148,0,0)#degminsec2decimal_degrees(148,15,0)
    145 
    146 #Main Domain of Hobart: first run JS 18/9/06
    147 d0 = [south, west]
    148 d1 = [south, east]
    149 d2 = [north, east]
    150 d3 = [north, west]
    151 
    152 polyAll, zone = convert_points_from_latlon_to_utm([d0, d1, d2, d3])
    153 refzone = zone
    154 #polyAll = [[500000, 5200000],[580000, 5200000],[580000, 5270000],[500000,5270000]] #073937
    155 #polyAll = [[500000, 5160000],[580000, 5160000],[580000, 5270000],[500000,5270000]] # try this is the above works
    156 #polyAll = [[510000, 5170000],[570000, 5170000],[570000, 5250000],[510000,5250000]] # try this is the above works
    157 #polyAll = [[520000, 5160000],[580000, 5170000],[580000, 5200000],[600000,5260000],[510000,5280000]] # try this is the above works
    158 polyAll = [[520000, 5170000],[580000, 5170000],[580000, 5200000],[590000,5240000],[520000,5260000]] # try this is the above works
    159 # bottom bright, topr, top, left
    160 
    161 # region to export for Alex to make bathymetry map
    162 e_min_area = 560000#500000 #490000
    163 e_max_area = 570000#560000 #580000
    164 n_min_area = 5200000#5240000#5160000
    165 n_max_area = 5230000#5260000#5270000
     169###############################
     170# Interior region definitions
     171###############################
    166172
    167173#Interior region - Hobart city area + Glenorchy, Kingston
    … …  
    188194
    189195# Bruny Island
    190 from anuga.utilities.polygon import read_polygon
    191196poly_bruny = read_polygon(polygondir+'bruny.csv')

  • anuga_work/production/hobart_2006/run_hobart.py

    r3669 r3671  
    6161
    6262# filenames
    63 onshore_dem_name = project.onshore_dem_name
    64 onshore_dem_name_25 = project.onshore_dem_name_25
     63onshore_offshore_dem_name = project.onshore_offshore_dem_name
    6564meshname = project.meshname+'.msh'
    6665source_dir = project.boundarydir
    … …  
    6867copied_files = False
    6968
    70 # creates DEM from asc data - 12.5m
    71 #convert_dem_from_ascii2netcdf(onshore_dem_name, use_cache=True, verbose=True)
    72 
    73 #creates pts file for onshore DEM - 12.5
    74 #dem2pts(onshore_dem_name,
    75 #        easting_min=project.eastingmin,
    76 #        easting_max=project.eastingmax,
    77 #        northing_min=project.northingmin,
    78 #        northing_max= project.northingmax,
    79 #        use_cache=True, verbose=True)
    80 
    8169# create DEM from asc data - 25m data
    82 convert_dem_from_ascii2netcdf(onshore_dem_name_25, use_cache=True, verbose=True)
    83 
    84 #creates pts file for onshore DEM - 25
    85 dem2pts(onshore_dem_name_25,
    86         easting_min=project.eastingmin25,
    87         easting_max=project.eastingmax25,
    88         northing_min=project.northingmin25,
    89         northing_max= project.northingmax25,
    90         use_cache=True, verbose=True)
    91 
    92 #combine_rectangular_points_files(project.onshore_dem_name + '.pts',
    93 #                                 project.onshore_dem_name_25 + '.pts',
    94 #                                 project.all_onshore_dem_name + '.pts')
    95 
    96 print 'adding data sets'
    97 G = Geospatial_data(file_name = project.offshore_dem_name_local1 + '.xya')+\
    98     Geospatial_data(file_name = project.offshore_dem_name_local2 + '.xya')+\
    99     Geospatial_data(file_name = project.offshore_dem_name_local3 + '.xya')+\
    100     Geospatial_data(file_name = project.offshore_dem_name_local4 + '.xya')+\
    101     Geospatial_data(file_name = project.offshore_dem_name_aho1 + '.xya')+\
    102     Geospatial_data(file_name = project.offshore_dem_name_aho2 + '.xya')+\
    103     Geospatial_data(file_name = project.offshore_dem_name_aho3 + '.xya')+\
    104     Geospatial_data(file_name = project.offshore_dem_name_aho4 + '.xya')+\
    105     Geospatial_data(file_name = project.offshore_dem_name_aho5 + '.xya')+\
    106     Geospatial_data(file_name = project.offshore_dem_name_aho6 + '.xya')+\
    107     Geospatial_data(file_name = project.offshore_dem_name_aho7 + '.xya')+\
    108     Geospatial_data(file_name = project.offshore_dem_name_aho8 + '.xya')+\
    109     Geospatial_data(file_name = project.offshore_dem_name_aho9 + '.xya')+\
    110     Geospatial_data(file_name = project.offshore_dem_name_aho10 + '.xya')+\
    111     Geospatial_data(file_name = project.offshore_dem_name_aho11 + '.xya')+\
    112     Geospatial_data(file_name = project.offshore_dem_name_aho12 + '.xya')+\
    113     Geospatial_data(file_name = project.offshore_dem_name_aho13 + '.xya')+\
    114     Geospatial_data(file_name = project.offshore_dem_name_aho14 + '.xya')+\
    115     Geospatial_data(file_name = project.offshore_dem_name_aho15 + '.xya')+\
    116     Geospatial_data(file_name = project.offshore_dem_name_aho16 + '.xya')+\
    117     Geospatial_data(file_name = project.onshore_dem_name_25 + '.pts')
     70convert_dem_from_ascii2netcdf(onshore_offshore_dem_name, use_cache=True, verbose=True)
     71
     72#creates pts file for combined DEM - 25
     73dem2pts(onshore_offshore_dem_name, use_cache=True, verbose=True)
     74
     75# create geospatial data set and export
     76G = Geospatial_data(file_name = project.onshore_offshore_dem_name + '.pts')
    11877G.export_points_file(project.combined_dem_name + '.pts')
    119 #Geospatial_data(file_name = project.all_onshore_dem_name + '.pts')
     78
    12079#----------------------------------------------------------------------------
    12180# Create the triangular mesh based on overall clipping polygon with a tagged
    … …  
    228187# 7 min square wave starting at 1 min, 6m high
    229188Bw = Time_boundary(domain = domain,
    230                    f=lambda t: [(60<t<480)*6, 0, 0])
     189                   f=lambda t: [(60<t<480)*10, 0, 0])
    231190
    232191# for MOST BC
    … …  
    236195# for testing
    237196domain.set_boundary( {'topr': Bd, 'left': Bd, 'top': Bd,
    238                       'bottom': Bd, 'bright': Bw} )
     197                      'bottom': Bw, 'bright': Bd} )
    239198
    240199#-------------------------------------------------------------------------------                                 
Note: See TracChangeset for help on using the changeset viewer.