Changeset 5785


Ignore:
Timestamp:
Sep 25, 2008, 1:34:14 PM (15 years ago)
Author:
kristy
Message:

Correlating with Janes updates

Location:
anuga_work/production/perth
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • anuga_work/production/perth/build_perth_250m.py

    r5782 r5785  
    11"""
    2 Script for building the elevation data to run tsunami inundation scenario
     2Script for building the elevation data to run a tsunami inundation scenario
    33for Perth, WA, Australia.
    44
    … …  
    77The run_perth.py is reliant on the output of this script.
    88
    9 Ole Nielsen and Duncan Gray, GA - 2005 and Jane Sexton, Nick Bartzis, GA - 2006
    109"""
    1110
    … …  
    2928
    3029# Application specific imports
    31 import project_250m   # Definition of file names and polygons
     30import project   # Definition of file names and polygons
    3231
    3332#------------------------------------------------------------------------------
    … …  
    5049# Fine pts file to be clipped to area of interest
    5150#-------------------------------------------------------------------------------
    52 print"project_250m.poly_all",project_250m.poly_all
    53 print"project_250m.combined_dir_name",project_250m.combined_dir_name
     51print "project_250m.poly_all", project_250m.poly_all
     52print "project_250m.combined_dir_name", project_250m.combined_dir_name
    5453
    5554# input elevation directory filenames
    … …  
    6362convert_dem_from_ascii2netcdf(onshore_in_dir_name, basename_out=onshore_dir_name, use_cache=True, verbose=True)
    6463
    65 #creates pts file for onshore DEM
     64# creates pts file for onshore DEM
    6665print "creates pts file for onshore DEM"
    6766dem2pts(onshore_dir_name ,use_cache=True,verbose=True)
    6867
    69 #-------------------------------------------------------------------------------
    70 # Combine datasets into project_250m.combined_dir_name
    71 #-------------------------------------------------------------------------------
    72 
     68# create onshore pts files
    7369print'create Geospatial data1 objects from topographies'
    7470G = Geospatial_data(file_name = onshore_dir_name + '.pts')
     71
     72#-------------------------------------------------------------------------------
     73# Combine, clip and export dataset
     74#-------------------------------------------------------------------------------
    7575
    7676print'clip combined geospatial object by bounding polygon'
    … …  
    8383#G_clipped.export_points_file(project_250m.combined_dir_name + '.txt') #Use for comparision in ARC
    8484
    85 import sys
    86 sys.exit()

  • anuga_work/production/perth/export_results_all.py

    r5781 r5785  
    1111"""
    1212
    13 import project, os
     13import project_250m, os
    1414import sys
    1515from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf, dem2pts
    … …  
    1818from os import sep
    1919
    20 directory = project.output_dir
     20directory = project_250m.output_dir
    2121
    2222#time_dir = '20080526_104946_run_final_0.6_test_kvanputt'
    2323#time_dir = '20080530_170833_run_final_0.6_exmouth_kvanputt'
    2424#time_dir1 = '20080815_103442_run_final_0.0_polyline_alpha0.1_kvanputt'
    25 time_dir1 = '20080912_154439_run_final_0.6_27255_alpha0.1_kvanputt'
     25time_dir = '20080924_123601_run_final_0_27283_250m_none_kvanputt'
    2626
    2727#cellsize = 20
    2828cellsize = 30
    29 #timestep = 0
     29timestep = 0
    3030#area = ['Geordie', 'Sorrento', 'Fremantle', 'Rockingham']
    31 area = 'All'
    32 which_area = area
    33 var = [1,2] # Absolute momentum and depth
     31#area = ['All']
     32which_area = 'all'
     33#var = [1,2] # Absolute momentum and depth
    3434#var = [2] # depth
    35 #var = [0,4] #stage and elevation
     35var = [4] #stage and elevation
    3636
    37 time_dirs = [time_dir1]
    38 for time_dir in time_dirs:
     37name1 = directory+time_dir+sep+project_250m.scenario_name
     38#name2 = directory+time_dir+sep+'sww2'+sep+project.scenario_name+'_time_39900_0' #need to get assistance on how to make this into anything
    3939
    40     name1 = directory+time_dir+sep+project.scenario_name
    41     #name2 = directory+time_dir+sep+'sww2'+sep+project.scenario_name+'_time_39900_0' #need to get assistance on how to make this into anything
    42 
    43     names = [name1]
    44     for name in names:
     40names = [name1]
     41for name in names:
    4542
    4643##        for which_area in area:
    … …  
    7168##                northing_max = project.ymaxRockingham
    7269
    73              
    74             for which_var in var:
    75                 if which_var == 0:  # Stage
    76                     outname = name + which_area + '_stage'
    77                     quantityname = 'stage'
     70         
     71        for which_var in var:
     72            if which_var == 0:  # Stage
     73                outname = name + which_area + '_stage'
     74                quantityname = 'stage'
    7875
    79                 if which_var == 1:  # Absolute Momentum
    80                     outname = name + which_area + '_momentum'
    81                     quantityname = '(xmomentum**2 + ymomentum**2)**0.5' 
     76            if which_var == 1:  # Absolute Momentum
     77                outname = name + which_area + '_momentum'
     78                quantityname = '(xmomentum**2 + ymomentum**2)**0.5' 
    8279
    83                 if which_var == 2:  # Depth
    84                     outname = name + which_area + '_depth'
    85                     quantityname = 'stage-elevation' 
     80            if which_var == 2:  # Depth
     81                outname = name + which_area + '_depth'
     82                quantityname = 'stage-elevation' 
    8683
    87                 if which_var == 3:  # Speed
    88                     outname = name + which_area + '_speed'
    89                     #quantityname = '(xmomentum**2 + ymomentum**2)**0.5/(stage-elevation+1.e-6/(stage-elevation))'  #Speed
    90                     quantityname = '(xmomentum**2 + ymomentum**2)**0.5/(stage-elevation+1.e-6)'  #Speed
    91                    
    92                 if which_var == 4:  # Elevation
    93                     outname = name + which_area + '_elevation'
    94                     quantityname = 'elevation'  #Elevation
     84            if which_var == 3:  # Speed
     85                outname = name + which_area + '_speed'
     86                #quantityname = '(xmomentum**2 + ymomentum**2)**0.5/(stage-elevation+1.e-6/(stage-elevation))'  #Speed
     87                quantityname = '(xmomentum**2 + ymomentum**2)**0.5/(stage-elevation+1.e-6)'  #Speed
     88               
     89            if which_var == 4:  # Elevation
     90                outname = name + which_area + '_elevation'
     91                quantityname = 'elevation'  #Elevation
    9592
    96                 else:
    97                     print 'start sww2dem',area
    98                     sww2dem(name, basename_out = outname,
    99                                 quantity = quantityname,
    100                                 #timestep = timestep,
    101                                 cellsize = cellsize,     
    102                                 reduction = max,
    103                                 verbose = True,
    104                                 format = 'asc')
     93            else:
     94                print 'start sww2dem',which_area
     95                sww2dem(name, basename_out = outname,
     96                            quantity = quantityname,
     97                            timestep = timestep,
     98                            cellsize = cellsize,     
     99                            reduction = max,
     100                            verbose = True,
     101                            format = 'asc')
    105102##
    106103##                else:

  • anuga_work/production/perth/project_250m.py

    r5782 r5785  
    11"""Common filenames and locations for elevation, meshes and outputs.
    2 This script is the heart of all scripts in folder
     2This script is the heart of all scripts in the folder
    33"""
    44#------------------------------------------------------------------------------
    … …  
    1818# Directory setup
    1919#------------------------------------------------------------------------------
    20 ##Note: INUNDATIONHOME is the inundation directory, not the data directory.
     20# Note: INUNDATIONHOME is the inundation directory, not the data directory.
    2121
    2222home = getenv('INUNDATIONHOME') + sep +'data'+sep #Sandpit's parent diruser = get_user_name()
    … …  
    2525host = get_host_name()
    2626
    27 ##time stuff
    28 time = strftime('%Y%m%d_%H%M%S',localtime()) #gets time for new dir
    29 gtime = strftime('%Y%m%d_%H%M%S',gmtime()) #gets time for new dir
     27# determines time for setting up output directories
     28time = strftime('%Y%m%d_%H%M%S',localtime())
     29gtime = strftime('%Y%m%d_%H%M%S',gmtime())
    3030build_time = time+'_build'
    3131run_time = time+'_run'
    32 print 'gtime: ', gtime
    3332
    3433#------------------------------------------------------------------------------
    … …  
    3635#------------------------------------------------------------------------------
    3736
    38 ##Changed to reflect the modeled community (also sets up the directory system)
     37# this section needs to be updated to reflect the modelled community.
     38# Note, the user needs to set up the directory system accordingly
    3939state = 'western_australia'
    4040scenario_name = 'perth'
    4141scenario = 'perth_tsunami_scenario'
    42 ##One or all can be changed each time the run_scenario script is excuted
    43 tide = 0 #0.6
    44 #event_number = 27255
     42
     43# Model specific parameters. One or all can be changed each time the
     44# run_scenario script is executed
     45tide = 0                #0.6
     46#event_number = 27255   # linked to hazard map
    4547event_number = 27283
    46 alpha = 0.1
    47 friction=0.01
    48 starttime=0
    49 finaltime=80000
     48alpha = 0.1             # smoothing parameter for mesh
     49friction=0.01           # manning's friction coefficient
     50starttime=0             
     51finaltime=80000         # final time for simulation
    5052
    5153interior_mesh = 'all' # Can have 'all' or 'none' for Phase 2 study
    5254
    53 setup='final'  ## Can replace with trial or basic, this action will thin the mesh
    54                ## so that the run script will take less time (hence less acurate)
     55setup='final'  # Final can be replaced with trial or basic.
     56               # Either will result in a coarser mesh that will allow a
     57               # faster, but less accurate, simulation.
    5558
    5659if setup =='trial':
    … …  
    7376# Output Filename
    7477#------------------------------------------------------------------------------
    75 ##Important to distiquish each run
     78# Important to distinguish each run - ensure str(user) is included!
     79# Note, the user is free to include as many parameters as desired
    7680dir_comment='_'+setup+'_'+str(tide)+'_'+str(event_number)+'_250m_' + str(interior_mesh) +'_'+str(user)
    7781
    7882#------------------------------------------------------------------------------
    79 # INPUT DATA
    80 #------------------------------------------------------------------------------
    81 
    82 ##ELEVATION DATA## - used in build_perth.py
    83 ## onshore data: format ascii grid with accompaning projection file
    84 grid_250m_2005 = 'grid'
    85 
    86 ##GAUGES## - used in get_timeseries.py
     83# Input Data
     84#------------------------------------------------------------------------------
     85
     86# elevation data used in build_perth.py
     87# onshore data: format ascii grid with accompanying projection file
     88onshore_name = 'grid'
     89
     90# gauges - used in get_timeseries.py
    8791gauge_name = 'perth.csv'
    8892gauge_name2 = 'thinned_MGA50.csv'
    8993
    90 ##BOUNDING POLYGON## -used in build_boundary.py and run_perth.py respectively
    91 #NOTE: when files are put together must be in sequence - for ease go clockwise!
    92 #Check the run_perth.py for boundary_tags
    93 ##thinned ordering file from Hazard Map: format index,latitude,longitude (with title)
     94# BOUNDING POLYGON - used in build_boundary.py and run_perth.py respectively
     95# NOTE: when files are put together the points must be in sequence - for ease go clockwise!
     96# Check the run_perth.py for boundary_tags
     97# thinned ordering file from Hazard Map: format is index,latitude,longitude (with title)
    9498order_filename = 'thinned_boundary_ordering.txt'
    95 ##landward bounding points
     99#landward bounding points
    96100landward = 'landward_bounding_polygon.txt'
    97101
    98102#------------------------------------------------------------------------------
    99 # OUTPUT ELEVATION DATA
    100 #------------------------------------------------------------------------------
    101 ##Output filename for elevation
    102 ## its a combination of all the data put together (utilisied in build_boundary)
     103# Output Elevation Data
     104#------------------------------------------------------------------------------
     105# Output filename for elevation
     106# this is a combination of all the data (utilisied in build_boundary)
    103107combined_name ='perth_combined_elevation_250m'
    104108
    … …  
    117121
    118122#------------------------------------------------------------------------------
    119 # Location of input and output Data
    120 #------------------------------------------------------------------------------
    121 ##where the input data sits
    122 onshore_in_dir_name = topographies_in_dir + grid_250m_2005
    123 
    124 ##where the output data sits
    125 onshore_dir_name = topographies_dir + grid_250m_2005
    126 
    127 ##where the combined file sits
     123# Location of input and output data
     124#------------------------------------------------------------------------------
     125# where the input data sits
     126onshore_in_dir_name = topographies_in_dir + onshore_name
     127
     128# where the output data sits
     129onshore_dir_name = topographies_dir + onshore_name
     130
     131# where the combined elevation file sits
    128132combined_dir_name = topographies_dir + combined_name
    129133
    130 ##where the mesh sits (this is created during the run_perth.py)
    131 meshes_dir_name = meshes_dir + scenario_name + interior_mesh + '.msh'
    132 
    133 #where the boundary order files sit (this is used within build_boundary.py)
     134# where the mesh sits (this is created during the run_perth.py)
     135meshes_dir_name = meshes_dir + scenario_name+ interior_mesh +'.msh'
     136
     137# where the boundary ordering files sit (this is used within build_boundary.py)
    134138order_filename_dir = boundaries_dir + order_filename
    135139
    136 #where the landward points of boundary extent sit (this is used within run_perth.py)
     140# where the landward points of boundary extent sit (this is used within run_perth.py)
    137141landward_dir = boundaries_dir + landward
    138142
    139 #where the event sts files sits (this is created during the build_boundary.py)
     143# where the event sts files sits (this is created during the build_boundary.py)
    140144boundaries_dir_event = boundaries_dir + str(event_number) + sep
    141145boundaries_dir_mux = muxhome
    142146
    143 #where the directory of the output filename sits
    144 output_build_time_dir = output_dir+build_time+dir_comment+sep #used for build_perth.py
    145 output_run_time_dir = output_dir+run_time+dir_comment+sep #used for run_perth.py
     147# where the directory of the output filename sits
     148output_build_time_dir = output_dir+build_time+dir_comment+sep   #used for build_perth.py
     149output_run_time_dir = output_dir+run_time+dir_comment+sep       #used for run_perth.py
    146150output_run_time_dir_name = output_run_time_dir + scenario_name  #Used by post processing
    147151
    148 #where the directory of the gauges sit
    149 gauges_dir_name = gauges_dir + gauge_name #used for get_timeseries.py
    150 gauges_dir_name2 = gauges_dir + gauge_name2 #used for get_timeseries.py
     152#w here the directory of the gauges sit
     153gauges_dir_name = gauges_dir + gauge_name       #used for get_timeseries.py
     154gauges_dir_name2 = gauges_dir + gauge_name2     #used for get_timeseries.py
    151155
    152156#------------------------------------------------------------------------------
    … …  
    154158#------------------------------------------------------------------------------
    155159
    156 #Land, to set the stage/water to be offcoast only
     160#Land, to set the initial stage/water to be offcoast only
    157161poly_mainland = read_polygon(polygons_dir+'initial_condition.csv')
    158162
    159 #Initial bounding polygon for data clipping
     163# Initial bounding polygon for data clipping
    160164poly_all = read_polygon(polygons_dir+'poly_all.csv')
    161165res_poly_all = 100000*res_factor
    162166
    163 #Area of Interest 1 (Fremantle)
     167# Area of Interest 1 (Fremantle)
    164168poly_aoi1 = read_polygon(polygons_dir+'CBD_coastal.csv')
    165169res_aoi1 = 500*res_factor
    166170
    167 #Area of Interest 2 (Rockingham)
     171# Area of Interest 2 (Rockingham)
    168172poly_aoi2 = read_polygon(polygons_dir+'rockingham_penguin.csv')
    169173res_aoi2 = 500*res_factor
    170174
    171 #Area of Interest 2 (garden Island)
     175# Area of Interest 2 (garden Island)
    172176poly_aoi2a = read_polygon(polygons_dir+'garden.csv')
    173177res_aoi2a= 500*res_factor
    174178
    175 #Area of Interest 3 (geordie bay - record of tsunami impact)
     179# Area of Interest 3 (geordie bay - record of tsunami impact)
    176180poly_aoi3 = read_polygon(polygons_dir+'geordie_bay.csv')
    177181res_aoi3 = 500*res_factor
    178182
    179 #Area of Interest 4 (sorrento - record of tsunami impact)
     183# Area of Interest 4 (sorrento - record of tsunami impact)
    180184poly_aoi4 = read_polygon(polygons_dir+'sorrento_gauge.csv')
    181185res_aoi4 = 500*res_factor
    182186
    183 #Area of Significance 1 (Garden Island and sand bank infront of Rockingham)
     187# Area of Significance 1 (Garden Island and sand bank infront of Rockingham)
    184188poly_aos1 = read_polygon(polygons_dir+'garden_rockingham.csv')
    185189res_aos1 = 1000*res_factor
    186190
    187 #Area of Significance 2 (incorporate coastline of rottnest)
     191# Area of Significance 2 (incorporate coastline of rottnest)
    188192poly_aos2 = read_polygon(polygons_dir+'rottnest_external.csv')
    189193res_aos2 = 1000*res_factor
    190194
    191 #Refined areas
    192 #Polygon designed to incorporate Dredge Area from Fremantle to
    193 #Rockingham the steep incline was making the mesh go to 0
     195# Refined areas
     196# Polygon designed to incorporate dredged area from Fremantle to
     197# Rockingham as the steep incline was making the elevation go to 0
    194198poly_aos3 = read_polygon(polygons_dir+'DredgeArea.csv')
    195199res_aos3 = 1000*res_factor
    196200
    197 #Shallow water 1
     201# Shallow water 1
    198202poly_sw1 = read_polygon(polygons_dir+'internal_h20mORd3km.csv')
    199203res_sw1 = 25000*res_factor
    200204
    201 #Deep water (land of Rottnest, ANUGA does not do donuts!)
     205# Deep water (land of Rottnest, ANUGA does not do donuts!)
    202206poly_dw1 = read_polygon(polygons_dir+'rottnest_internal.csv')
    203207res_dw1 = 100000*res_factor
    … …  
    215219    print'Mesh = none'
    216220    interior_regions = []
    217    
     221     
    218222trigs_min = number_mesh_triangles(interior_regions, poly_all, res_poly_all)
    219223print 'min estimated number of triangles', trigs_min
    220224   
    221 
    222225#------------------------------------------------------------------------------
    223226# Clipping regions for export to asc and regions for clipping data
     227# Final inundation maps should only be created in regions of the finest mesh
    224228#------------------------------------------------------------------------------
    225229

  • anuga_work/production/perth/run_perth_250m.py

    r5782 r5785  
    1 """Script for running tsunami inundation scenario for Perth, WA, Australia.
     1"""Script for running a tsunami inundation scenario for Perth, WA, Australia.
    22
    33The scenario is defined by a triangular mesh created from project_250m.polygon,
    4 the elevation data is combiled into a pts file through build_perth.py
    5 and a simulated tsunami is generated through an sts file from build_boundary.py
    6 
    7 Input: sts file (build_boundary.py for resepective event)
     4the elevation data is compiled into a pts file through build_perth.py
     5and a simulated tsunami is generated through an sts file from build_boundary.py.
     6
     7Input: sts file (build_boundary.py for respective event)
    88       pts file (build_perth.py)
    9        information from project_250m file
     9       information from project file
    1010Outputs: sww file stored in project_250m.output_run_time_dir
    1111The export_results_all.py and get_timeseries.py is reliant
    … …  
    4848   
    4949# Application specific imports
    50 import project_250m  # Definition of file names and polygons
     50import project  # Definition of file names and polygons
    5151numprocs = 1
    5252myid = 0
    … …  
    9090    N = len(urs_bounding_polygon)-1
    9191
    92     #boundary tags refer to project_250m.landward 4 points equals 5 segments start at N
     92    # boundary tags refer to project_250m.landward 4 points equals 5 segments start at N
    9393    boundary_tags={'back': [N+1,N+2,N+3], 'side': [N,N+4], 'ocean': range(N)}
    9494
    95    
    9695    #--------------------------------------------------------------------------
    9796    # Create the triangular mesh based on overall clipping polygon with a tagged
    … …  
    10099    #--------------------------------------------------------------------------
    101100
    102     #IMPORTANT don't cache create_mesh_from_region and Domain(mesh....) as it
     101    # IMPORTANT don't cache create_mesh_from_region and Domain(mesh....) as it
    103102    # causes problems with the ability to cache set quantity which takes alot of times
    104103       
    … …  
    132131    print 'Setup initial conditions'
    133132
    134     #following sets the stage/water to be offcoast only
     133    # sets the initial stage in the offcoast region only
    135134    IC = Polygon_function( [(project_250m.poly_mainland, 0)], default = kwargs['tide'],
    136135                             geo_reference = domain.geo_reference)
    … …  
    161160    domain.set_name(kwargs['scenario_name'])
    162161    domain.set_datadir(kwargs['output_dir'])
    163     domain.set_default_order(2) # Apply second order scheme
    164     domain.set_minimum_storable_height(0.01) # Don't store anything less than 1cm
     162    domain.set_default_order(2)                 # Apply second order scheme
     163    domain.set_minimum_storable_height(0.01)    # Don't store anything less than 1cm
    165164    domain.set_store_vertices_uniquely(False)
    166165    domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
    … …  
    206205        domain.write_boundary_statistics(tags = 'ocean')
    207206
     207    # these outputs should be checked with the resultant inundation map
    208208    x, y = domain.get_maximum_inundation_location()
    209209    q = domain.get_maximum_inundation_elevation()
    210 
    211210    print 'Maximum runup observed at (%.2f, %.2f) with elevation %.2f' %(x,y,q)
    212211
    213     print 'That took %.2f seconds' %(time.time()-t0)
     212    print 'Simulation took %.2f seconds' %(time.time()-t0)
    214213
    215214    #kwargs 'completed' must be added to write the final parameters to file
    … …  
    232231    kwargs['alpha'] = project_250m.alpha
    233232    kwargs['friction']=project_250m.friction
    234        
    235 
     233     
    236234    run_model(**kwargs)
    237235     
Note: See TracChangeset for help on using the changeset viewer.