Changeset 3271


Ignore:
Timestamp:
Jul 4, 2006, 11:04:04 AM (19 years ago)
Author:
sexton
Message:

update pt hedland scripts to write output results to perlite 2 - note, must change INUNDATIONHOME path variable to \perlite\cit\2\cit\data

Location:
production/pt_hedland_2006
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • production/pt_hedland_2006/get_building_inundation.py

    r2896 r3271  
    2929
    3030# Inputs
     31#timestampdir = '' # MSL
    3132#timestampdir = '' # HAT
    3233timestampdir = '' # LAT

  • production/pt_hedland_2006/project.py

    r3190 r3271  
    1010from coordinate_transforms.redfearn import degminsec2decimal_degrees
    1111from time import localtime, strftime
    12 
    1312from geospatial_data import *
    1413
    15 
    16 
    17 
    1814#Making assumptions about the location of scenario data
     15state = 'western_australia'
    1916scenario_dir_name = 'pt_hedland_tsunami_scenario_2006'
    2017
    2118# onshore data from 30m DTED level 2
    22 onshore_name = 'pt_hedland_onshore_30m_dted' # get from Neil/Ingo (DEM or topo data)
     19onshore_name_dted = 'pt_hedland_onshore_30m_dted' # get from Neil/Ingo (DEM or topo data)
     20onshore_name_dli = 'pt_hedland_onshore_20m_dli' # get from Neil/Ingo (DEM or topo data)
     21
    2322# offshore data from GA digitised charts
    2423offshore_name1 = 'pt_hedland_offshore_points'
     24
    2525# offshore data from AHO fairsheets
    2626offshore_name2 = 'pt_hedland_offshore_points_fairsheet'
     27
    2728# coastline developed from 30m DTED
    2829coast_name = 'pt_hedland_coastline_points.xya'
    … …  
    3637
    3738if sys.platform == 'win32':
    38     home = getenv('INUNDATIONHOME')
     39    home = getenv('INUNDATIONHOME') #Sandpit's parent dir
    3940#    python_home = getenv('PWD')     
    40 #    home = environ['INUNDATIONHOME']     #Sandpit's parent dir
    4141    #user = basename(getenv('USERPROFILE'))
    4242    #print 'USER:', user
    43 else:   
    44     home = getenv('INUNDATIONHOME', sep+'d'+sep+'cit'+sep+'1'+sep+'cit'+sep+'risk_assessment_methods_project'+sep+'inundation')     
     43else:
     44    # original
     45    #home = getenv('INUNDATIONHOME', sep+'d'+sep+'cit'+sep+'1'+sep+'cit'+sep+'risk_assessment_methods_project'+sep+'inundation')
     46    # update to perlite 2
     47    home = getenv('INUNDATIONHOME', sep+'d'+sep+'cit'+sep+'2'+sep+'cit'+sep+'inundation'+sep+'data')     
    4548    user = getenv('LOGNAME')
    4649    print 'USER:', user
    … …  
    4851#Derive subdirectories and filenames
    4952time = strftime('%Y%m%d_%H%M%S',localtime()) #gets time for new dir
    50 #print 'home', home
    51 outputtimedir = home+sep+scenario_dir_name+sep+'output'+sep+time+sep
     53outputtimedir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'outputs'+sep+time+sep
     54
    5255#print 'outputtimedir', outputtimedir
    53 meshdir = home+sep+scenario_dir_name+sep+'meshes'+sep
    54 datadir = home+sep+scenario_dir_name+sep+'topographies'+sep
    55 gaugedir = home+sep+scenario_dir_name+sep+'gauges'+sep
    56 polygondir = home+sep+scenario_dir_name+sep+'polygons'+sep
    57 boundarydir = home+sep+scenario_dir_name+sep+'boundaries'+sep
    58 #output dir without time
    59 outputdir = home+sep+scenario_dir_name+sep+'output'+sep
    60 tidedir = home+sep+scenario_dir_name+sep+'tide_data'+sep
     56#meshdir = home+sep+scenario_dir_name+sep+'meshes'+sep
     57#datadir = home+sep+scenario_dir_name+sep+'topographies'+sep
     58#gaugedir = home+sep+scenario_dir_name+sep+'gauges'+sep
     59#polygondir = home+sep+scenario_dir_name+sep+'polygons'+sep
     60#boundarydir = home+sep+scenario_dir_name+sep+'boundaries'+sep
     61#outputdir = home+sep+scenario_dir_name+sep+'outputs'+sep
     62#tidedir = home+sep+scenario_dir_name+sep+'tide_data'+sep
    6163
    62 #print'bound', boundarydir
     64meshdir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'meshes'+sep
     65datadir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'topographies'+sep
     66gaugedir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'gauges'+sep
     67polygondir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'polygons'+sep
     68boundarydir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'boundaries'+sep
     69outputdir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'outputs'+sep
     70tidedir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'tide_data'+sep
    6371
    64 #gauge_filename = gaugedir + 'onslow_gauges.xya'
    65 #for MOST
    66 #gauge_filename = gaugedir + 'pt_hedland_gauges.xya'
    6772gauge_filename = gaugedir + 'gauge_location_port_hedland.csv'
    6873buildings_filename = gaugedir + 'pt_hedland_res.csv'
    … …  
    7277coast_filename = datadir + coast_name
    7378
    74 # boundary source data
    75 #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
    76 
    77 #print 'name', __name__
    78 #print 'path', __file__
    79 #codedir = getcwd()+sep
    80 
    81 #project_code_name = __name__
    82                            
    83 #project_code_dir_name = __file__
    84 
    8579meshname = meshdir + basename
    86 
    87 onshore_dem_name = datadir + onshore_name
    88 
     80onshore_dem_name = datadir + onshore_name_dted
    8981offshore_dem_name1 = datadir + offshore_name1
    9082offshore_dem_name2 = datadir + offshore_name2
    91 
    9283combined_dem_name = datadir + 'pt_hedland_combined_elevation'
    93 
    9484outputname = outputtimedir + basename  #Used by post processing
    9585
    … …  
    10797
    10898# region to export (used from export_results.py)
    109 
    11099e_min_area = 633000
    111100e_max_area = 690000
    … …  
    113102n_max_area = 7761000
    114103
    115 refzone = 50 # confirm with Hamish
     104refzone = 50
    116105
    117106# bounding polygon provided by Hamish
    … …  
    155144
    156145poly_region = [j0, j1, j2, j3]
    157 
    158 coast_buffer_file = datadir+'pts2ascii_test.xya'
    159 G = Geospatial_data(file_name=coast_buffer_file,delimiter=' ')
    160 poly_coast = list(G.get_data_points())
    161 #print 'get_data_points()',G.get_data_points()
    162 #print 'get_',poly_region
    163 
    164 
    165 

  • production/pt_hedland_2006/run_pt_hedland.py

    r3251 r3271  
    1 """Script for running a tsunami inundation scenario for Port Hedland, WA, Australia.
     1"""Script for running a tsunami inundation scenario for Onslow, WA, Australia.
    22
    33Source data such as elevation and boundary data is assumed to be available in
    … …  
    55The output sww file is stored in project.outputtimedir
    66
    7 The scenario is defined by a triangular mesh created from project.polygon and
    8 the elevation data.
    9 
    10 Ole Nielsen and Duncan Gray, GA - 2005, Nick Bartzis and Jane Sexton, GA - 2006
     7The scenario is defined by a triangular mesh created from project.polygon,
     8the elevation data and a simulated submarine landslide.
     9
     10Ole Nielsen and Duncan Gray, GA - 2005 and Nick Bartzis, GA - 2006
    1111"""
    12 
    13 
    14 #------------------------------------------------------------------------------
    15 # Import necessary modules
    16 #------------------------------------------------------------------------------
     12#-------------------------------------------------------------------------------# Import necessary modules
     13#-------------------------------------------------------------------------------
    1714
    1815# Standard modules
    19 import os
    20 
    2116from os import sep
    2217from os.path import dirname, basename
    23 
    2418import time
    2519
    … …  
    3024from pyvolution.combine_pts import combine_rectangular_points_files
    3125from pyvolution.pmesh2domain import pmesh_to_domain_instance
    32 #from geospatial_data import add_points_files
    33 
    34 # Application specific imports
    35 import project                 # Definition of file names and polygons
    36 from smf import slump_tsunami  # Function for submarine mudslide
    37 
    3826from shutil import copy
    3927from os import mkdir, access, F_OK
    40 
    4128from geospatial_data import *
    4229import sys
    4330from pyvolution.util import Screen_Catcher
    4431
    45 #------------------------------------------------------------------------------
     32# Application specific imports
     33import project                 # Definition of file names and polygons
     34
     35
     36#-------------------------------------------------------------------------------
    4637# Preparation of topographic data
    4738#
    … …  
    4940# Do for coarse and fine data
    5041# Fine pts file to be clipped to area of interest
    51 #------------------------------------------------------------------------------
     42#-------------------------------------------------------------------------------
    5243
    5344# filenames
    … …  
    5849source_dir = project.boundarydir
    5950
    60 
    61 
    62 
    63 #import sys; sys.exit()
    64 
    65 # creates copy of code in output dir if dir doesn't exist
     51## creates copy of code in output dir if dir doesn't exist
    6652if access(project.outputtimedir,F_OK) == 0 :
    6753    mkdir (project.outputtimedir)
    … …  
    6955copy (__file__, project.outputtimedir + basename(__file__))
    7056
    71 
    72 #normal screen output is stored in
     57print 'project.outputtimedir',project.outputtimedir
     58
     59##normal screen output is stored in
    7360screen_output_name = project.outputtimedir + "screen_output.txt"
    7461screen_error_name = project.outputtimedir + "screen_error.txt"
    7562
    7663#used to catch screen output to file
    77 #sys.stdout = Screen_Catcher(screen_output_name)
    78 #sys.stderr = Screen_Catcher(screen_error_name)
    79 
    80 
    81 '''
     64sys.stdout = Screen_Catcher(screen_output_name)
     65sys.stderr = Screen_Catcher(screen_error_name)
     66
    8267# fine data (clipping the points file to smaller area)
    8368# creates DEM from asc data
    … …  
    9075        northing_min=project.northingmin,
    9176        northing_max= project.northingmax,
    92         use_cache=True, 
     77        use_cache=True,
    9378        verbose=True)
    9479
    95 print'create G1'
     80print 'create G1'
    9681G1 = Geospatial_data(file_name = project.offshore_dem_name1 + '.xya')
    9782
    98 print'create G2'
     83print 'create G2'
    9984G2 = Geospatial_data(file_name = project.offshore_dem_name2 + '.xya')
    10085
    101 print'create G3'
     86print 'create G3'
    10287G3 = Geospatial_data(file_name = project.onshore_dem_name + '.pts')
    10388
    104 print'add G1+G2+G3'
     89print 'add G1+G2+G3'
    10590G = G1 + G2 + G3
    10691
    107 print'export G'
     92print 'export G'
    10893G.export_points_file(project.combined_dem_name + '.pts')
    10994
    110 '''
    111 #------------------------------------------------------------------------------
     95
     96#-------------------------------------------------------------------------------                                 
    11297# Create the triangular mesh based on overall clipping polygon with a tagged
    11398# boundary and interior regions defined in project.py along with
    11499# resolutions (maximal area of per triangle) for each polygon
    115 #------------------------------------------------------------------------------
     100#-------------------------------------------------------------------------------
    116101
    117102from pmesh.mesh_interface import create_mesh_from_regions
    118103
    119104region_res = 100000
    120 coast_res = 2500
    121 pt_hedland_res = 500
     105coast_res = 25000
     106pt_hedland_res = 5000
    122107# derive poly_coast from project.coast_name using alpha_shape
     108#interior_regions = [[project.poly_pt_hedland, pt_hedland_res],
    123109interior_regions = [[project.poly_pt_hedland, pt_hedland_res],
    124                     [project.poly_coast, coast_res]]
    125 #                    [project.poly_region, region_res]]
     110                    [project.poly_region, region_res]]
    126111
    127112print 'number of interior regions', len(interior_regions)
    … …  
    139124        count += 1
    140125
    141 figname = 'pt_hedland_polys'
    142 print'pt_hedland_polys'
    143 #plot_polygons([project.polyAll, project.poly_pt_hedland, project.poly_region],
    144 plot_polygons([project.poly_coast],
     126if sys.platform == 'win32':
     127    #figname = project.outputtimedir + 'pt_hedland_polys'
     128    figname = 'pt_hedland_polys'
     129    plot_polygons([project.polyAll,project.poly_pt_hedland,project.poly_region],
    145130              figname,
    146               verbose = True)
     131              verbose = True)   
    147132
    148133if count == 0:
    … …  
    152137    print 'check %s in production directory' %figname
    153138    import sys; sys.exit()
    154 
     139   
    155140print 'start create mesh from regions'
    156141from caching import cache
    157142_ = cache(create_mesh_from_regions,
    158143          project.polyAll,
     144#          project.poly_region,
    159145#          {'boundary_tags': {'right': [0], 'bottomright': [1],
    160146#                             'bottomleft': [2], 'left': [3], 'top': [4]},
    161           {'boundary_tags': {'topright': [0], 'top': [1],'topleft': [2], 'left': [3],
    162                              'bottomleft': [4], 'bottomright': [5], 'right': [6]},
    163            'maximum_triangle_area': 10000000,
     147#          {'boundary_tags': {'topright': [0], 'top': [1],'topleft': [2], 'left': [3],
     148#                             'bottomleft': [4], 'bottomright': [5], 'right': [6]},
     149          {'boundary_tags': {'topright': [0],'topleft': [1], 'left': [2],
     150                             'bottomleft': [3], 'bottomright': [4], 'right': [5]},
     151           'maximum_triangle_area': 500000,
    164152           'filename': meshname,           
    165153           'interior_regions': interior_regions},
    166           verbose = True)
    167 
    168 
    169 #------------------------------------------------------------------------------
     154          verbose = True, evaluate=True)
     155
     156
     157#-------------------------------------------------------------------------------                                 
    170158# Setup computational domain
    171 #------------------------------------------------------------------------------
    172 
    173 '''
    174 domain = pmesh_to_domain_instance(meshname, Domain,
    175                                   use_cache = False,
    176                                   verbose = True)
    177 '''
    178 
     159#-------------------------------------------------------------------------------                                 
    179160domain = Domain(meshname, use_cache = False, verbose = True)
    180161
     162print domain.statistics()
    181163print 'Number of triangles = ', len(domain)
    182164print 'The extent is ', domain.get_extent()
    … …  
    187169domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
    188170
    189 #------------------------------------------------------------------------------
     171#-------------------------------------------------------------------------------                                 
    190172# Setup initial conditions
    191 #------------------------------------------------------------------------------
     173#-------------------------------------------------------------------------------
    192174
    193175tide = 0.
     176#high
     177#tide = 3.6
     178#low
     179#tide = -3.9
    194180
    195181domain.set_quantity('stage', tide)
    … …  
    204190                    )
    205191
    206 #print 'have sent quantities OK - now exiting'
    207 #import sys; sys.exit()
    208 
    209 #------------------------------------------------------------------------------
     192#-------------------------------------------------------------------------------                                 
    210193# Setup boundary conditions (all reflective)
    211 #------------------------------------------------------------------------------
     194#-------------------------------------------------------------------------------
    212195print 'start ferret2sww'
    213196# skipped as results in file SU-AU_clipped is correct for all WA
    … …  
    235218       'fail_on_NaN': False,
    236219       'inverted_bathymetry': True},
    237       #evaluate = True,
    238        verbose = True,
    239       dependencies= source_dir + project.boundary_basename + '.sww')
     220       evaluate = True,
     221       verbose = True,)
    240222
    241223
    … …  
    247229Bd = Dirichlet_boundary([tide,0,0])
    248230
    249 # 7 min square wave starting at 1 min, 6m high
    250 Bw = Time_boundary(domain = domain,
    251                    f=lambda t: [(60<t<480)*6, 0, 0])
    252 
    253 #domain.set_boundary( {'top': Bf, 'left': Br, 'bottomleft': Br, 'bottomright': Br, 'right': Br} )
    254 #
    255 domain.set_boundary( {'topright': Bf, 'top': Bf,'topleft': Bf, 'left':  Br,
    256                              'bottomleft': Br, 'bottomright': Br, 'right': Br})
     231domain.set_boundary( {'topright': Bf,'topleft': Bf, 'left':  Bd,
     232                             'bottomleft': Bd, 'bottomright': Bd, 'right': Bd})
    257233                             
    258 #------------------------------------------------------------------------------
     234#-------------------------------------------------------------------------------                                 
    259235# Evolve system through time
    260 #------------------------------------------------------------------------------
     236#-------------------------------------------------------------------------------
    261237import time
    262238t0 = time.time()
    263239
    264 for t in domain.evolve(yieldstep = 240, finaltime = 7200):
    265     domain.write_time()
    266     domain.write_boundary_statistics(tags = 'top')     
    267 
    268 for t in domain.evolve(yieldstep = 120, finaltime = 12600
    269                        ,skip_initial_step = True):
    270     domain.write_time()
    271     domain.write_boundary_statistics(tags = 'top')     
    272 
    273 for t in domain.evolve(yieldstep = 60, finaltime = 19800
    274                        ,skip_initial_step = True):
    275     domain.write_time()
    276     domain.write_boundary_statistics(tags = 'top')     
     240for t in domain.evolve(yieldstep = 240, finaltime = 12240):
     241    domain.write_time()
     242    domain.write_boundary_statistics(tags = 'topleft')     
     243
     244for t in domain.evolve(yieldstep = 120, finaltime = 15600
     245                       ,skip_initial_step = True):
     246    domain.write_time()
     247    domain.write_boundary_statistics(tags = 'topleft')     
     248
     249for t in domain.evolve(yieldstep = 60, finaltime = 22020
     250                       ,skip_initial_step = True):
     251    domain.write_time()
     252    domain.write_boundary_statistics(tags = 'topleft')     
    277253   
    278 for t in domain.evolve(yieldstep = 120, finaltime = 25200
    279                        ,skip_initial_step = True):
    280     domain.write_time()
    281     domain.write_boundary_statistics(tags = 'top')     
     254for t in domain.evolve(yieldstep = 120, finaltime = 27060
     255                       ,skip_initial_step = True):
     256    domain.write_time()
     257    domain.write_boundary_statistics(tags = 'topleft')     
    282258
    283259for t in domain.evolve(yieldstep = 240, finaltime = 36000
    284260                       ,skip_initial_step = True):
    285261    domain.write_time()
    286     domain.write_boundary_statistics(tags = 'top')     
     262    domain.write_boundary_statistics(tags = 'topleft')     
    287263 
    288264print 'That took %.2f seconds' %(time.time()-t0)
Note: See TracChangeset for help on using the changeset viewer.