Changeset 6382


Ignore:
Timestamp:
Feb 23, 2009, 4:16:51 PM (14 years ago)
Author:
kristy
Message:

working new script

Location:
anuga_work/production/australia_ph2/broome
Files:
3 deleted
4 edited

Legend:

Unmodified
Added
Removed

  • anuga_work/production/australia_ph2/broome/build_elevation.py

    r6339 r6382  
    1 """Build the elevation data to run a tsunami scenario
    2 for Broome, WA, Australia for input into the Ph2 Hazard Map.
     1"""Build the elevation data to run a tsunami inundation scenario
     2for Busselton, WA, Australia.
    33
    44Input: elevation data from project.py
    … …  
    2424
    2525# Application specific imports
    26 import project   # Definition of file names and polygons
     26from setup_model import project   # Definition of file names and polygons
    2727
    2828
    … …  
    5050# Create Geospatial data from ASCII files
    5151geospatial_data = {}
    52 ##for filename in project.ascii_grid_filenames:
    53 ##    absolute_filename = join(project.topographies_folder, filename)
    54 ##    convert_dem_from_ascii2netcdf(absolute_filename,
    55 ##                                  basename_out=absolute_filename,
    56 ##                                  use_cache=True,
    57 ##                                  verbose=True)
    58 ##    dem2pts(absolute_filename, use_cache=True, verbose=True)
    59 ##
    60 ##    geospatial_data[filename] = Geospatial_data(file_name=absolute_filename+
    61 ##                                                '.pts', verbose=True)
     52if not project.ascii_grid_filenames == []:
     53    for filename in project.ascii_grid_filenames:
     54        absolute_filename = join(project.topographies_folder, filename)
     55        convert_dem_from_ascii2netcdf(absolute_filename,
     56                                      basename_out=absolute_filename,
     57                                      use_cache=True,
     58                                      verbose=True)
     59        dem2pts(absolute_filename, use_cache=True, verbose=True)
     60
     61        geospatial_data[filename] = Geospatial_data(file_name=absolute_filename+'.pts',
     62                                                    verbose=True)
    6263
    6364# Create Geospatial data from TXT files
    64 for filename in project.point_filenames:
    65     absolute_filename = join(project.topographies_folder, filename)
    66     geospatial_data[filename] = Geospatial_data(file_name=absolute_filename,
    67                                                 verbose=True)
    68 
     65if not project.point_filenames == []:
     66    for filename in project.point_filenames:
     67        absolute_filename = join(project.topographies_folder, filename)
     68        geospatial_data[filename] = Geospatial_data(file_name=absolute_filename,
     69                                                    verbose=True)
    6970
    7071#-------------------------------------------------------------------------------

  • anuga_work/production/australia_ph2/broome/file_length.py

    r6339 r6382  
    77Returns: number of lines in file
    88"""
     9
    910def file_length(in_file):
     11    '''Function to return the number of lines in a file.
     12
     13    in_file: Path to the file to get number of lines in.
     14   
     15    Returns: number of lines in file
     16    '''
     17   
    1018    fid = open(in_file)
    1119    data = fid.readlines()

  • anuga_work/production/australia_ph2/broome/project.py

    r6346 r6382  
    1 """Common filenames and locations for elevation, meshes and outputs.
    2 This script is the heart of all scripts in the folder
    31"""
    4 #------------------------------------------------------------------------------
    5 # Import necessary modules
    6 #------------------------------------------------------------------------------
     2This file contains all your file and directory definitions
     3for elevation, meshes and outputs.
     4"""
    75
    86import os
    9 import sys
    10 from os.path import join
    11 from os import sep, getenv
     7from anuga.utilities.system_tools import get_user_name, get_host_name
    128from time import localtime, strftime, gmtime
    13 from anuga.utilities.polygon import read_polygon, number_mesh_triangles
    14 from anuga.utilities.system_tools import get_user_name, get_host_name
    15 
    16 #------------------------------------------------------------------------------
     9from os.path import join, exists
     10
     11
     12#-------------------------------------------------------------------------------
    1713# Directory setup
    18 #------------------------------------------------------------------------------
    19 # Note: INUNDATIONHOME is the inundation directory, not the data directory.
    20 
    21 home = getenv('INUNDATIONHOME')+sep+'data'+sep # Absolute path for data folder
    22 muxhome = getenv('MUXHOME')
    23 user = get_user_name()
    24 host = get_host_name()
    25 
    26 # determines time for setting up output directories
    27 time = strftime('%Y%m%d_%H%M%S',localtime())
    28 gtime = strftime('%Y%m%d_%H%M%S',gmtime())
    29 build_time = time+'_build'
    30 run_time = time+'_run'
     14#-------------------------------------------------------------------------------
    3115
    3216# this section needs to be updated to reflect the modelled community.
    … …  
    3418state = 'australia_ph2'
    3519scenario_name = 'broome'
    36 
    37 #------------------------------------------------------------------------------
     20scenario_folder = scenario_name
     21
     22#-------------------------------------------------------------------------------
    3823# Initial Conditions
    39 #------------------------------------------------------------------------------
    40 # Model specific parameters. One or all can be changed each time the
    41 # run_scenario script is executed
    42 tide = 0
    43 event_number = 27255 # Java 9.3 worst case for Perth
     24#-------------------------------------------------------------------------------
     25
     26# Model specific parameters.
     27# One or all can be changed each time the run_model script is executed
     28tide = 0                # difference between MSL and HAT
     29event_number = 27255    # the event number or the mux file name
    4430alpha = 0.1             # smoothing parameter for mesh
    45 friction = 0.01           # manning's friction coefficient
    46 starttime = 0             
    47 finaltime = 1000         # final time for simulation
    48 
    49 setup = 'trial'  # Final can be replaced with trial or basic.
    50                # Either will result in a coarser mesh that will allow a
    51                # faster, but less accurate, simulation.
    52 
    53 if setup =='trial':
    54     print'trial'
    55     scale_factor=100
    56     time_thinning=96
    57     yieldstep=240
    58 if setup =='basic':
    59     print'basic'
    60     scale_factor=4
    61     time_thinning=12
    62     yieldstep=120
    63 if setup =='final':
    64     print'final'
    65     scale_factor=1
    66     time_thinning=4
    67     yieldstep=60
    68 
    69 
    70 #------------------------------------------------------------------------------
    71 # Output Filename
    72 #------------------------------------------------------------------------------
    73 # Important to distinguish each run - ensure str(user) is included!
    74 # Note, the user is free to include as many parameters as desired
    75 output_comment= ('_' + setup + '_' + str(tide)+ '_' + str(event_number) +
    76                  '_' + str(user))
    77 
    78 #------------------------------------------------------------------------------
     31friction=0.01           # manning's friction coefficient
     32starttime=0             # start time for simulation
     33finaltime=1000          # final time for simulation
     34
     35setup = 'trial'         # This can be one of three values
     36                        #    trial - coarsest mesh, fast
     37                        #    basic - coarse mesh
     38                        #    final - fine mesh, slowest
     39
     40#-------------------------------------------------------------------------------
     41# Output filename
     42#
     43# Your output filename should be unique between different runs on different data.
     44# The list of items below will be used to create a file in your output directory.
     45# Your user name and time+date will be automatically added.  For example,
     46#     [setup, tide, event_number]
     47# will result in a filename like
     48#     20090212_091046_run_final_0_27283_rwilson
     49#-------------------------------------------------------------------------------
     50
     51output_comment = [setup, tide, event_number]
     52
     53#-------------------------------------------------------------------------------
    7954# Input Data
    80 #------------------------------------------------------------------------------
     55#-------------------------------------------------------------------------------
     56
    8157# ELEVATION DATA
    8258# Used in build_elevation.py
    83 
    84 ### Format for ascii grids, as produced in ArcGIS + a projection file
    85 ##ascii_grid_filenames = ['grid250m'] # 250m grid 2005
     59# Format for ascii grids, as produced in ArcGIS + a projection file
     60ascii_grid_filenames = [] # Grid data
    8661
    8762# Format for point is x,y,elevation (with header)
    88 point_filenames = ['grid250m.txt'] # 250m grid 2005
     63point_filenames = ['grid.txt']            # 250m grid 2005
    8964
    9065# BOUNDING POLYGON - for data clipping and estimate of triangles in mesh
    … …  
    9267# Format for points easting,northing (no header)
    9368bounding_polygon_filename = 'bounding_polygon.csv'
     69bounding_polygon_maxarea = 100000
     70
     71# INTERIOR REGIONS -  for designing the mesh
     72# Used in run_model.py
     73# Format for points easting,northing (no header)                   
     74interior_regions_data = []
     75
     76# LAND - used to set the initial stage/water to be offcoast only
     77# Used in run_model.py.  Format for points easting,northing (no header)
     78land_initial_conditions_filename = []
    9479
    9580# GAUGES - for creating timeseries at a specific point
    96 # Used in get_timeseries.py
     81# Used in get_timeseries.py. 
    9782# Format easting,northing,name,elevation (with header)
    98 ##gauges_filename = 'gauges.csv'
    99 
    100 # BOUNDING POLYGON
    101 # used in build_boundary.py and run_model.py respectively
     83gauges_filename = 'gauges.csv'
     84
     85# BUILDINGS EXPOSURE - for identifying inundated houses
     86# Used in run_building_inundation.py
     87# Format latitude,longitude etc (geographic)
     88building_exposure_filename = 'busselton_res_clip.csv' # from NEXIS
     89
     90# BOUNDING POLYGON - used in build_boundary.py and run_model.py respectively
    10291# NOTE: when files are put together the points must be in sequence
    10392# For ease go clockwise!
    … …  
    112101landward_boundary_filename = 'landward_boundary.csv'
    113102
    114 #------------------------------------------------------------------------------
     103# MUX input filename.
     104# If a meta-file from EventSelection is used, set 'multi-mux' to True.
     105# If a single MUX stem filename (*.grd) is used, set 'multi-mux' to False.
     106##mux_input_filename = event_number # to be found in event_folder
     107                                    # (ie boundaries/event_number/)
     108##multi_mux = False
     109mux_input_filename = 'event.list'
     110multi_mux = True
     111
     112
     113################################################################################
     114################################################################################
     115####         NOTE: NOTHING WOULD NORMALLY CHANGE BELOW THIS POINT.          ####
     116################################################################################
     117################################################################################
     118
     119# Get system user and host names.
     120# These values can be used to distinguish between two similar runs by two
     121# different users or runs by the same user on two different machines.
     122user = get_user_name()
     123host = get_host_name()
     124
     125# Environment variable names.
     126# The inundation directory, not the data directory.
     127ENV_INUNDATIONHOME = 'INUNDATIONHOME'
     128
     129# Path to MUX data
     130ENV_MUXHOME = 'MUXHOME'
     131
     132#-------------------------------------------------------------------------------
    115133# Output Elevation Data
    116 #------------------------------------------------------------------------------
     134#-------------------------------------------------------------------------------
     135
    117136# Output filename for elevation
    118137# this is a combination of all the data generated in build_elevation.py
    119 combined_elevation_basename = scenario_name + '_combined_elevation_new'
    120 
    121 #------------------------------------------------------------------------------
     138combined_elevation_basename = scenario_name + '_combined_elevation'
     139
     140#-------------------------------------------------------------------------------
    122141# Directory Structure
    123 #------------------------------------------------------------------------------
    124 anuga_folder = join(home, state, scenario_name, 'anuga')
     142#-------------------------------------------------------------------------------
     143
     144# determines time for setting up output directories
     145time = strftime('%Y%m%d_%H%M%S', localtime())
     146gtime = strftime('%Y%m%d_%H%M%S', gmtime())
     147build_time = time + '_build'
     148run_time = time + '_run_'
     149
     150# create paths generated from environment variables.
     151home = join(os.getenv(ENV_INUNDATIONHOME), 'data') # Absolute path for data folder
     152muxhome = os.getenv(ENV_MUXHOME)
     153   
     154# check various directories/files that must exist
     155anuga_folder = join(home, state, scenario_folder, 'anuga')
    125156topographies_folder = join(anuga_folder, 'topographies')
    126157polygons_folder = join(anuga_folder, 'polygons')
    127158boundaries_folder = join(anuga_folder, 'boundaries')
    128159output_folder = join(anuga_folder, 'outputs')
    129 gauges_folder = join(anuga_folder,'gauges')
     160gauges_folder = join(anuga_folder, 'gauges')
    130161meshes_folder = join(anuga_folder, 'meshes')
    131 
    132 #------------------------------------------------------------------------------
     162event_folder = join(boundaries_folder, str(event_number))
     163
     164# MUX data files
     165# Directory containing the MUX data files to be used with EventSelection.
     166mux_data_folder = join(muxhome, 'mux')
     167
     168#-------------------------------------------------------------------------------
    133169# Location of input and output data
    134 #------------------------------------------------------------------------------
     170#-------------------------------------------------------------------------------
     171
     172# Convert the user output_comment to a string for run_model.py
     173output_comment = ('_'.join([str(x) for x in output_comment if x != user])
     174                  + '_' + user)
    135175
    136176# The absolute pathname of the all elevation, generated in build_elevation.py
    … …  
    138178
    139179# The absolute pathname of the mesh, generated in run_model.py
    140 meshes = join(meshes_folder, scenario_name) + 'new.msh'
    141 
    142 # The absolute pathname for the urs order points, used within build_boundary.py
     180meshes = join(meshes_folder, scenario_name) + '.msh'
     181
     182# The pathname for the urs order points, used within build_urs_boundary.py
    143183urs_order = join(boundaries_folder, urs_order_filename)
    144184
    … …  
    147187landward_boundary = join(boundaries_folder, landward_boundary_filename)
    148188
    149 event_folder = join(boundaries_folder, str(event_number))
    150                    
    151189# The absolute pathname for the .sts file, generated in build_boundary.py
    152190event_sts = join(event_folder, scenario_name)
    … …  
    162200# The absolute pathname for the gauges file
    163201# Used for get_timeseries.py
    164 ##gauges = join(gauges_folder, gauges_filename)       
    165 
    166 #------------------------------------------------------------------------------
    167 # Reading polygons and creating interior regions
    168 #------------------------------------------------------------------------------
    169 
    170 # Initial bounding polygon for data clipping
    171 bounding_polygon = read_polygon(join(polygons_folder,
    172                                      bounding_polygon_filename))
    173 bounding_maxarea = 100000*scale_factor
    174 
    175 interior_regions = []
    176 
    177 # Estimate the number of triangles                     
    178 trigs_min = number_mesh_triangles(interior_regions,
    179                                   bounding_polygon, bounding_maxarea)
    180 print 'min estimated number of triangles', trigs_min
    181    
    182 
     202gauges = join(gauges_folder, gauges_filename)       
     203
     204# The absolute pathname for the building file
     205# Used for run_building_inundation.py
     206building_exposure = join(gauges_folder, building_exposure_filename)
     207
     208# full path to where MUX files (or meta-files) live
     209mux_input = join(event_folder, mux_input_filename)
     210

  • anuga_work/production/australia_ph2/broome/run_model.py

    r6340 r6382  
    2222# Standard modules
    2323import os
     24import os.path
    2425import time
     26from time import localtime, strftime, gmtime
    2527
    2628# Related major packages
     29from Scientific.IO.NetCDF import NetCDFFile
     30import Numeric as num
     31
    2732from anuga.interface import create_domain_from_regions
    2833from anuga.interface import Transmissive_stage_zero_momentum_boundary
    … …  
    3641from anuga.shallow_water.data_manager import start_screen_catcher
    3742from anuga.shallow_water.data_manager import copy_code_files
     43from anuga.shallow_water.data_manager import urs2sts
    3844from anuga.utilities.polygon import read_polygon, Polygon_function
    39    
     45
    4046# Application specific imports
    41 import project  # Definition of file names and polygons
     47from setup_model import project
     48import build_urs_boundary as bub
    4249
    43 
    44 #------------------------------------------------------------------------------
     50#-------------------------------------------------------------------------------
    4551# Copy scripts to time stamped output directory and capture screen
    4652# output to file. Copy script must be before screen_catcher
    47 #------------------------------------------------------------------------------
     53#-------------------------------------------------------------------------------
     54
    4855copy_code_files(project.output_run, __file__,
    49                 os.path.dirname(project.__file__)+os.sep+\
    50                 project.__name__+'.py' )
     56                os.path.join(os.path.dirname(project.__file__),
     57                             project.__name__+'.py'))
    5158start_screen_catcher(project.output_run, 0, 1)
    5259
    53 
    54 #------------------------------------------------------------------------------
     60#-------------------------------------------------------------------------------
    5561# Create the computational domain based on overall clipping polygon with
    5662# a tagged boundary and interior regions defined in project.py along with
    5763# resolutions (maximal area of per triangle) for each polygon
    58 #------------------------------------------------------------------------------
     64#-------------------------------------------------------------------------------
     65
    5966print 'Create computational domain'
     67
     68# Create the STS file
     69print 'project.mux_data_folder=%s' % project.mux_data_folder
     70if not os.path.exists(project.event_sts + '.sts'):
     71    bub.build_urs_boundary(project.mux_input_filename, project.event_sts)
    6072
    6173# Read in boundary from ordered sts file
    … …  
    7082
    7183# Number of boundary segments
    72 N = len(event_sts)-1
     84num_ocean_segments = len(event_sts) - 1
     85# Number of landward_boundary points
     86num_land_points = file_length(project.landward_boundary)
    7387
    74 # Number of landward_boundary points
    75 M = file_length(project.landward_boundary)
    76                
    7788# Boundary tags refer to project.landward_boundary
    7889# 4 points equals 5 segments start at N
    79 boundary_tags={'back': range(N+1, N+M),
    80                'side': [N, N+M],
    81                'ocean': range(N)}
     90boundary_tags={'back': range(num_ocean_segments+1,
     91                             num_ocean_segments+num_land_points),
     92               'side': [num_ocean_segments,
     93                        num_ocean_segments+num_land_points],
     94               'ocean': range(num_ocean_segments)}
    8295
    8396# Build mesh and domain
    … …  
    95108domain.set_minimum_storable_height(0.01)    # Don't store depth less than 1cm
    96109
     110#-------------------------------------------------------------------------------
     111# Setup initial conditions
     112#-------------------------------------------------------------------------------
    97113
    98 #------------------------------------------------------------------------------
    99 # Setup initial conditions
    100 #------------------------------------------------------------------------------
    101114print 'Setup initial conditions'
    102115
    103116# Set the initial stage in the offcoast region only
    104 ##IC = Polygon_function(project.land_initial_conditions,
    105 ##                      default=project.tide,
    106 ##                      geo_reference=domain.geo_reference)
    107 domain.set_quantity('stage', 0, use_cache=True, verbose=True)
     117if project.land_initial_conditions:
     118    IC = Polygon_function(project.land_initial_conditions,
     119                          default=project.tide,
     120                          geo_reference=domain.geo_reference)
     121else:
     122    IC = 0
     123domain.set_quantity('stage', IC, use_cache=True, verbose=True)
    108124domain.set_quantity('friction', project.friction)
    109125domain.set_quantity('elevation',
    … …  
    113129                    alpha=project.alpha)
    114130
     131#-------------------------------------------------------------------------------
     132# Setup boundary conditions
     133#-------------------------------------------------------------------------------
    115134
    116 #------------------------------------------------------------------------------
    117 # Setup boundary conditions
    118 #------------------------------------------------------------------------------
    119135print 'Set boundary - available tags:', domain.get_boundary_tags()
    120136
    121137Br = Reflective_boundary(domain)
    122138Bt = Transmissive_stage_zero_momentum_boundary(domain)
     139Bd = Dirichlet_boundary([project.tide, 0, 0])
    123140Bf = Field_boundary(project.event_sts+'.sts',
    124141                    domain, mean_stage=project.tide,
    … …  
    130147
    131148domain.set_boundary({'back': Br,
    132                      'side': Bt,
     149                     'side': Bd,
    133150                     'ocean': Bf})
    134151
     152#-------------------------------------------------------------------------------
     153# Evolve system through time
     154#-------------------------------------------------------------------------------
    135155
    136 #------------------------------------------------------------------------------
    137 # Evolve system through time
    138 #------------------------------------------------------------------------------
    139156t0 = time.time()
    140157for t in domain.evolve(yieldstep=project.yieldstep,
    … …  
    144161    print domain.boundary_statistics(tags='ocean')
    145162
    146 print 'Simulation took %.2f seconds' %(time.time()-t0)
     163print 'Simulation took %.2f seconds' % (time.time()-t0)
Note: See TracChangeset for help on using the changeset viewer.