Changeset 5790


Ignore:
Timestamp:
Sep 26, 2008, 11:25:31 AM (16 years ago)
Author:
kristy
Message:

Updated scripts to reflect perth. Also added the 250m scripts and used more than one polygon for initial condition

Location:
anuga_work/production/carnarvon
Files:
5 added
5 edited

Legend:

Unmodified
Added
Removed

  • anuga_work/production/carnarvon/build_carnarvon.py

    r5755 r5790  
    1 """Script for running tsunami inundation scenario for Dampier, WA, Australia.
     1"""
     2Script for building the elevation data to run a tsunami inundation scenario
     3for geraldton, WA, Australia.
    24
    3 Source data such as elevation and boundary data is assumed to be available in
    4 directories specified by project.py
    5 The output sww file is stored in project.output_time_dir
     5Input: elevation data from project.py
     6Output: pts file stored in project.topographies_dir
     7The run_geraldton.py is reliant on the output of this script.
    68
    7 The scenario is defined by a triangular mesh created from project.polygon,
    8 the elevation data and a simulated submarine landslide.
    9 
    10 Ole Nielsen and Duncan Gray, GA - 2005 and Jane Sexton, Nick Bartzis, GA - 2006
    119"""
    1210
    … …  
    2725from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf, dem2pts
    2826from anuga.geospatial_data.geospatial_data import *
    29 from anuga.shallow_water.data_manager import start_screen_catcher, copy_code_files,store_parameters
     27from anuga.shallow_water.data_manager import start_screen_catcher, copy_code_files
    3028
    3129# Application specific imports
    … …  
    4240start_screen_catcher(project.output_build_time_dir)
    4341
    44 print 'USER: ', project.user
     42print 'USER:    ', project.user
    4543
    4644#-------------------------------------------------------------------------------
    … …  
    5149# Fine pts file to be clipped to area of interest
    5250#-------------------------------------------------------------------------------
    53 print"project.poly_all",project.poly_all
    54 print"project.combined_dir_name",project.combined_dir_name
     51print "project.poly_all", project.poly_all
     52print "project.combined_dir_name", project.combined_dir_name
    5553
    56 # topography directory filenames
     54# input elevation directory filenames
    5755onshore_in_dir_name = project.onshore_in_dir_name
    5856coast_in_dir_name = project.coast_in_dir_name
    … …  
    6159offshore_in_dir_name2 = project.offshore_in_dir_name2
    6260
     61# output elevation directory filenames
    6362onshore_dir_name = project.onshore_dir_name
    6463coast_dir_name = project.coast_dir_name
    65 
    6664offshore_dir_name = project.offshore_dir_name
    6765offshore_dir_name1 = project.offshore_dir_name1
    … …  
    7068
    7169# creates DEM from asc data
    72 print "creates DEMs from ascii data"
     70print "creates DEMs from asc data"
    7371convert_dem_from_ascii2netcdf(onshore_in_dir_name, basename_out=onshore_dir_name, use_cache=True, verbose=True)
    7472
    75 #creates pts file for onshore DEM
     73# creates pts file for onshore DEM
    7674print "creates pts file for onshore DEM"
    77 dem2pts(onshore_dir_name,
    78 #        easting_min=project.eastingmin,
    79 #        easting_max=project.eastingmax,
    80 #        northing_min=project.northingmin,
    81 #        northing_max= project.northingmax,
    82         use_cache=True,
    83         verbose=True)
     75dem2pts(onshore_dir_name ,use_cache=True,verbose=True)
    8476
    85 #creates pts file for island DEM
    86 #dem2pts(island_dir_name, use_cache=True, verbose=True)
     77# create onshore pts files
     78print'create Geospatial data1 objects from topographies'
     79G1 = Geospatial_data(file_name = onshore_dir_name + '.pts')
    8780
    88 print'create Geospatial data1 objects from topographies'
    89 G1 = Geospatial_data(file_name = onshore_dir_name + '.pts',verbose=True)
    90 G2 = Geospatial_data(file_name = coast_in_dir_name + '.txt',verbose=True)
     81# create coastal and offshore txt files
     82print'create Geospatial data4 objects from topographies'
     83G_coast = Geospatial_data(file_name = coast_in_dir_name)
     84print'create Geospatial data5 objects from topographies'
     85G_off1 = Geospatial_data(file_name = offshore_in_dir_name)
     86print'create Geospatial data6 objects from topographies'
     87G_off2 = Geospatial_data(file_name = offshore_in_dir_name1)
     88print'create Geospatial data7 objects from topographies'
     89G_off3 = Geospatial_data(file_name = offshore_in_dir_name2)
    9190
    92 G_off = Geospatial_data(file_name = offshore_in_dir_name + '.txt',verbose=True)
    93 G_off1 = Geospatial_data(file_name = offshore_in_dir_name1 + '.txt',verbose=True)
    94 G_off2 = Geospatial_data(file_name = offshore_in_dir_name2 + '.txt',verbose=True)
     91#-------------------------------------------------------------------------------
     92# Combine, clip and export dataset
     93#-------------------------------------------------------------------------------
     94
     95print 'add all bathymetry files'
     96G_off = G_off1 + G_off2 + G_off3
     97print 'clipping data to coast'
     98G_ocean = G_off.clip(project.poly_ocean, verbose=True) #mainland
     99G_ocean_b = G_ocean.clip_outside(project.poly_island1, verbose=True) #bernier
     100G_ocean_bd = G_ocean_b.clip_outside(project.poly_island2, verbose=True) #dorne
     101
     102G1_clip_m = G1.clip(project.poly_mainland, verbose=True)
     103G1_clip_b = G1.clip(project.poly_island1, verbose=True)
     104G1_clip_d = G1.clip(project.poly_island2, verbose=True)
     105
    95106print'add all geospatial objects'
    96 G = G1 + G2 + G_off + G_off1 + G_off2
    97 
    98 print'clip combined geospatial object by bounding polygon'
    99 G_clipped = G.clip(project.poly_all)
     107G = G1_clip_m+ G1_clip_b+G1_clip_d+ G_coast + G_ocean_bd
    100108
    101109print'export combined DEM file'
    102110if access(project.topographies_dir,F_OK) == 0:
    103111    mkdir (project.topographies_dir)
    104 G_clipped.export_points_file(project.combined_dir_name + '.txt')
     112G.export_points_file(project.combined_dir_name + '.pts')
     113#G.export_points_file(project.combined_dir_name + '.txt') #Use for comparision in ARC
    105114
    106 print'project.combined_dir_name + .txt',project.combined_dir_name + '.txt'
    107 
    108 ###-------------------------------------------------------------------------
    109 ### Convert URS to SWW file for boundary conditions
    110 ###-------------------------------------------------------------------------
    111 ##print 'starting to create boundary conditions'
    112 ##from anuga.shallow_water.data_manager import urs2sww, urs_ungridded2sww
    113 ##
    114 ##boundaries_in_dir_name = project.boundaries_in_dir_name
    115 ##print 'boundaries_in_dir_name',project.boundaries_in_dir_name
    116 ##
    117 ##
    118 ###import sys; sys.exit()
    119 ##
    120 ##urs_ungridded2sww(project.boundaries_in_dir_name, project.boundaries_in_dir_name,
    121 ##                  verbose=True, mint=4000, maxt=80000, zscale=1)
    122 ##

  • anuga_work/production/carnarvon/export_results.py

    r5755 r5790  
    77
    88
    9 time_dir = '20080908_154507_run_final_0.0_polyline_alpha0.1_kvanputt'
     9#time_dir = '20080908_154507_run_final_0.0_polyline_alpha0.1_kvanputt'
     10time_dir = '20080909_123144_run_final_0.0_polyline_alpha0.1_kvanputt'
    1011
    1112

  • anuga_work/production/carnarvon/get_gauges.py

    r5755 r5790  
    1414        compress,zeros,fabs,take
    1515    fid = NetCDFFile(filename+'.sts', 'r')    #Open existing file for read
     16    permutation = fid.variables['permutation'][:]
    1617    x = fid.variables['x'][:]+fid.xllcorner   #x-coordinates of vertices
    1718    y = fid.variables['y'][:]+fid.yllcorner   #y-coordinates of vertices
    … …  
    1920    time=fid.variables['time'][:]+fid.starttime
    2021    elevation=fid.variables['elevation'][:]
    21     print 'hello', fid.variables.keys()
    22 
    23    
     22       
    2423    basename='sts_gauge'
    2524    quantity_names=['stage','xmomentum','ymomentum']
    … …  
    2928
    3029    maxname = 'max_sts_stage.csv'
    31     fid_max = open(project.boundaries_dir+sep+maxname,'w')
    32     s = 'x, y, max_stage \n'
     30    fid_max = open(project.boundaries_dir_event+sep+maxname,'w')
     31    s = 'index, x, y, max_stage \n'
    3332    fid_max.write(s)   
    3433    for j in range(len(x)):
    35         locx=int(x[j])
    36         locy=int(y[j])
     34        index= permutation[j]
    3735        stage = quantities['stage'][:,j]
    3836        xmomentum = quantities['xmomentum'][:,j]
    3937        ymomentum = quantities['ymomentum'][:,j]
    4038
    41         s = '%.6f, %.6f, %.6f\n' %(x[j], y[j], max(stage))
     39        s = '%d, %.6f, %.6f, %.6f\n' %(index, x[j], y[j], max(stage))
    4240        fid_max.write(s)
    4341       
    44         fid_sts = open(project.boundaries_dir+sep+basename+'_'+str(locx)+'_'+str(locy)+'.csv', 'w')
     42        fid_sts = open(project.boundaries_dir_event+sep+basename+'_'+ str(index)+'.csv', 'w')
    4543        s = 'time, stage, xmomentum, ymomentum \n'
    4644        fid_sts.write(s)
    … …  
    5654    return quantities,elevation,time
    5755
    58 quantities,elevation,time=get_sts_gauge_data(os.path.join(project.boundaries_dir,project.scenario_name),verbose=False)
     56quantities,elevation,time=get_sts_gauge_data(os.path.join(project.boundaries_dir_event,project.scenario_name),verbose=False)
    5957
    6058print len(elevation), len(quantities['stage'][0,:])

  • anuga_work/production/carnarvon/project.py

    r5755 r5790  
    1 # -*- coding: cp1252 -*-
    2 """Common filenames and locations for topographic data, meshes and outputs.
     1"""Common filenames and locations for elevation, meshes and outputs.
     2This script is the heart of all scripts in the folder
    33"""
     4#------------------------------------------------------------------------------
     5# Import necessary modules
     6#------------------------------------------------------------------------------
    47
    58from os import sep, environ, getenv, getcwd
    … …  
    811from time import localtime, strftime, gmtime
    912from anuga.utilities.polygon import read_polygon, plot_polygons, is_inside_polygon, number_mesh_triangles
    10 #from anuga.coordinate_transforms.redfearn import degminsec2decimal_degrees, convert_points_from_latlon_to_utm
    1113from anuga.utilities.system_tools import get_user_name, get_host_name
    1214from anuga.shallow_water.data_manager import urs2sts,create_sts_boundary
    1315from anuga.utilities.polygon import read_polygon, plot_polygons, polygon_area, is_inside_polygon
    1416
    15 # file and system info
    16 #---------------------------------
    17 #codename = 'project.py'
     17#------------------------------------------------------------------------------
     18# Directory setup
     19#------------------------------------------------------------------------------
     20# Note: INUNDATIONHOME is the inundation directory, not the data directory.
    1821
    1922home = getenv('INUNDATIONHOME') + sep +'data'+sep #Sandpit's parent diruser = get_user_name()
    … …  
    2225host = get_host_name()
    2326
    24 # INUNDATIONHOME is the inundation directory, not the data directory.
    25 
    26 #time stuff
    27 time = strftime('%Y%m%d_%H%M%S',localtime()) #gets time for new dir
    28 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())
    2930build_time = time+'_build'
    3031run_time = time+'_run'
    31 print 'gtime: ', gtime
    32 
    33 #Making assumptions about the location of scenario data
     32
     33#------------------------------------------------------------------------------
     34# Initial Conditions
     35#------------------------------------------------------------------------------
     36
     37# this section needs to be updated to reflect the modelled community.
     38# Note, the user needs to set up the directory system accordingly
    3439state = 'western_australia'
    3540scenario_name = 'carnarvon'
    3641scenario = 'carnarvon_tsunami_scenario'
    3742
    38 
    39 tide = 0.0 #1.0
    40 
    41 alpha = 0.1
    42 friction=0.01
    43 starttime=0
    44 finaltime=1000
    45 export_cellsize=25
    46 setup='final'
    47 source='polyline'
     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
     47event_number = 27283
     48alpha = 0.1             # smoothing parameter for mesh
     49friction=0.01           # manning's friction coefficient
     50starttime=0             
     51finaltime=80000         # final time for simulation
     52
     53setup='final'  # Final can be replaced with trial or basic.
     54               # Either will result in a coarser mesh that will allow a
     55               # faster, but less accurate, simulation.
    4856
    4957if setup =='trial':
    … …  
    6371    yieldstep=60
    6472
    65 dir_comment='_'+setup+'_'+str(tide)+'_'+str(source)+'_'+ 'alpha' +str(alpha)+'_'+str(user)
    66 
    67 # onshore data provided by WA DLI
    68 onshore_name = 'bathy250_clipland' # original
    69 
    70 # AHO + DPI data + colin French coastline
    71 coast_name = 'DPI_coastlineP'
    72 offshore_name = 'Shark_Bay_Clip'
    73 offshore_name1 = 'XYAHD_Clip'
    74 offshore_name2 = 'DPI_data'
    75 
    76 
    77 #final topo name
    78 combined_name = scenario_name+'_combined_elevation'
    79 combined_smaller_name = scenario_name+'_combined_elevation_smaller'
    80 
     73#------------------------------------------------------------------------------
     74# Output Filename
     75#------------------------------------------------------------------------------
     76# Important to distinguish each run - ensure str(user) is included!
     77# Note, the user is free to include as many parameters as desired
     78dir_comment='_'+setup+'_'+str(tide)+'_'+str(event_number)+'_'+ 'alpha' +str(alpha)+'_'+str(user)
     79
     80#------------------------------------------------------------------------------
     81# Input Data
     82#------------------------------------------------------------------------------
     83
     84# elevation data used in build_carnarvon.py
     85# onshore data: format ascii grid with accompanying projection file
     86onshore_name = 'dem_onshore'
     87# coastline: format x,y,elevation (with title)
     88coast_name = 'coastline.txt'
     89# bathymetry: format x,y,elevation (with title)
     90offshore_name = 'Shark_Bay_Clip.txt'
     91offshore_name1 = 'XYAHD_Clip.txt'
     92offshore_name2 = 'DPI_data.txt'
     93
     94# gauges - used in get_timeseries.py
     95gauge_name = 'carnarvon.csv'
     96gauge_name2 = 'thinned_MGA50.csv'
     97
     98# BOUNDING POLYGON - used in build_boundary.py and run_carnarvon.py respectively
     99# NOTE: when files are put together the points must be in sequence - for ease go clockwise!
     100# Check the run_carnarvon.py for boundary_tags
     101# thinned ordering file from Hazard Map: format is index,latitude,longitude (with title)
     102order_filename = 'thinned_boundary_ordering.txt'
     103#landward bounding points
     104landward = 'landward_bounding_polygon.txt'
     105
     106#------------------------------------------------------------------------------
     107# Output Elevation Data
     108#------------------------------------------------------------------------------
     109# Output filename for elevation
     110# this is a combination of all the data (utilisied in build_boundary)
     111combined_name ='carnarvon_combined_elevation'
     112combined_smaller_name = 'carnarvon_combined_elevation_smaller'
     113
     114#------------------------------------------------------------------------------
     115# Directory Structure
     116#------------------------------------------------------------------------------
    81117anuga_dir = home+state+sep+scenario+sep+'anuga'+sep
    82 
    83118topographies_in_dir = home+state+sep+scenario+sep+'elevation_final'+sep+'points'+sep
    84119topographies_dir = anuga_dir+'topographies'+sep
    85 
    86 # input topo file location
     120polygons_dir = anuga_dir+'polygons'+sep
     121tide_dir = anuga_dir+'tide_data'+sep
     122boundaries_dir = anuga_dir+'boundaries'+ sep
     123output_dir = anuga_dir+'outputs'+sep
     124gauges_dir = anuga_dir+'gauges'+sep
     125meshes_dir = anuga_dir+'meshes'+sep
     126
     127#------------------------------------------------------------------------------
     128# Location of input and output data
     129#------------------------------------------------------------------------------
     130# where the input data sits
    87131onshore_in_dir_name = topographies_in_dir + onshore_name
    88 
    89132coast_in_dir_name = topographies_in_dir + coast_name
    90133offshore_in_dir_name = topographies_in_dir + offshore_name
    … …  
    92135offshore_in_dir_name2 = topographies_in_dir + offshore_name2
    93136
     137# where the output data sits
    94138onshore_dir_name = topographies_dir + onshore_name
    95 
    96139coast_dir_name = topographies_dir + coast_name
    97140offshore_dir_name = topographies_dir + offshore_name
    … …  
    99142offshore_dir_name2 = topographies_dir + offshore_name2
    100143
    101 #final topo files
     144# where the combined elevation file sits
    102145combined_dir_name = topographies_dir + combined_name
    103 #combined_time_dir_name = topographies_time_dir + combined_name
    104146combined_smaller_name_dir = topographies_dir + combined_smaller_name
    105 #combined_time_dir_final_name = topographies_time_dir + combined_final_name
    106 
    107 meshes_dir = anuga_dir+'meshes'+sep
    108 meshes_dir_name = meshes_dir + scenario_name
    109 
    110 polygons_dir = anuga_dir+'polygons'+sep
    111 tide_dir = anuga_dir+'tide_data'+sep
    112 
    113 
    114 #boundaries_source = '1'
    115    
    116 if source=='polyline':
    117     boundaries_name = 'perth_3103_28052008' #polyline gun
    118     boundaries_in_dir = anuga_dir+'boundaries'+sep+'urs'+sep+'polyline'+sep+'1_10000'+sep
    119 
    120 if source=='test':
    121     boundaries_name = 'other' #polyline
    122     boundaries_in_dir = anuga_dir+'boundaries'+sep
    123 
    124 
    125 #boundaries locations
    126 boundaries_in_dir_name = boundaries_in_dir + boundaries_name
    127 boundaries_dir = anuga_dir+'boundaries'+sep
    128 boundaries_dir_name = boundaries_dir + scenario_name # what it creates???
     147
     148# where the mesh sits (this is created during the run_carnarvon.py)
     149meshes_dir_name = meshes_dir + scenario_name+'.msh'
     150
     151# where the boundary ordering files sit (this is used within build_boundary.py)
     152order_filename_dir = boundaries_dir + order_filename
     153
     154# where the landward points of boundary extent sit (this is used within run_carnarvon.py)
     155landward_dir = boundaries_dir + landward
     156
     157# where the event sts files sits (this is created during the build_boundary.py)
     158boundaries_dir_event = boundaries_dir + str(event_number) + sep
    129159boundaries_dir_mux = muxhome
    130160
    131 #output locations
    132 output_dir = anuga_dir+'outputs'+sep
    133 output_build_time_dir = anuga_dir+'outputs'+sep+build_time+dir_comment+sep
    134 output_run_time_dir = anuga_dir+'outputs'+sep+run_time+dir_comment+sep
     161# where the directory of the output filename sits
     162output_build_time_dir = output_dir+build_time+dir_comment+sep   #used for build_carnarvon.py
     163output_run_time_dir = output_dir+run_time+dir_comment+sep       #used for run_carnarvon.py
    135164output_run_time_dir_name = output_run_time_dir + scenario_name  #Used by post processing
    136165
    137 vertex_filename = output_run_time_dir + 'mesh_vertex.csv'
    138 
    139 #gauges
    140 gauge_name = 'carnarvon.csv'
    141 
    142 gauges_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'gauges'+sep
    143 beach_gauges = gauges_dir + 'beach_gauges.csv'
    144 gauges_dir_name = gauges_dir + gauge_name
    145 
    146 ##buildings_filename = gauges_dir + 'Perth_resA.csv'
    147 ##buildings_filename_out = 'Perth_res_Project_modified.csv'
    148 
    149 ###############################
     166#w here the directory of the gauges sit
     167gauges_dir_name = gauges_dir + gauge_name       #used for get_timeseries.py
     168gauges_dir_name2 = gauges_dir + gauge_name2     #used for get_timeseries.py
     169
     170#------------------------------------------------------------------------------
    150171# Interior region definitions
    151 ###############################
    152 
    153 #Initial bounding polygon for data clipping
     172#------------------------------------------------------------------------------
     173
     174#Land, to set the initial stage/water to be offcoast only
     175#Land and Ocean to clip data
     176poly_mainland=read_polygon(polygons_dir +'land_initial_condition.csv')
     177poly_island1=read_polygon(polygons_dir +'land_initial_condition_bernier.csv')
     178poly_island2=read_polygon(polygons_dir +'land_initial_condition_dorne.csv')
     179poly_ocean=read_polygon(polygons_dir +'ocean_initial_condition.csv')
     180
     181# Initial bounding polygon for data clipping
    154182poly_all = read_polygon(polygons_dir+'poly_all.csv')
    155183res_poly_all = 100000*res_factor
    156184
    157 #Polygon designed
    158 poly_internal_20 = read_polygon(polygons_dir+'Carnarvon20m.csv')
    159 res_internal_20 = 25000*res_factor
    160 
    161 #Polygon designed to
    162 poly_internal_5 = read_polygon(polygons_dir+'Carnarvon5m.csv')
    163 res_internal_5 = 500*res_factor
    164 
    165 #Polygon designed to
    166 poly_internal_10 = read_polygon(polygons_dir+'Carnarvon10m.csv')
    167 res_internal_10 = 1500*res_factor
    168 
    169 #Polygon designed to incorporate
    170 poly_island_20 = read_polygon(polygons_dir+'Island20m.csv')
    171 res_island_20 = 50000*res_factor
    172 
    173 
    174 interior_regions = [[poly_internal_20,res_internal_20],[poly_internal_5,res_internal_5]
    175                      ,[poly_island_20,res_island_20],[poly_internal_10,res_internal_10]]
     185# Area of Interest 1 (carnarvon)
     186poly_aoi1 = read_polygon(polygons_dir+'Carnarvon5m.csv')
     187res_aoi1 = 500*res_factor
     188
     189# Area of Significance 1 (carnarvon)
     190poly_aos1 = read_polygon(polygons_dir+'Carnarvon10m.csv')
     191res_aos1 = 1500*res_factor
     192
     193# Shallow water 1
     194poly_sw1 = read_polygon(polygons_dir+'Carnarvon20m.csv')
     195res_sw1 = 25000*res_factor
     196
     197# Shallow water 2
     198poly_sw2 = read_polygon(polygons_dir+'Island20m.csv')
     199res_sw2 = 50000*res_factor
     200
     201# Combined all regions, must check that all are included!
     202interior_regions = [[poly_aoi1,res_aoi1],[poly_aos1,res_aos1]
     203                     ,[poly_sw1,res_sw1],[poly_sw2,res_sw2]]
    176204
    177205   
    178206trigs_min = number_mesh_triangles(interior_regions, poly_all, res_poly_all)
    179 print 'min number triangles', trigs_min
     207print 'min estimated number of triangles', trigs_min
    180208   
    181 
    182 poly_mainland = read_polygon(polygons_dir+'initial_condition.csv')
    183 
    184 ###################################################################
     209#------------------------------------------------------------------------------
    185210# Clipping regions for export to asc and regions for clipping data
    186 ###################################################################
    187 
    188 
     211# Final inundation maps should only be created in regions of the finest mesh
     212#------------------------------------------------------------------------------
     213
     214# carnarvon CBD extract ascii grid
     215##xminCBD =
     216##xmaxCBD =
     217##yminCBD =
     218##ymaxCBD =

  • anuga_work/production/carnarvon/run_carnarvon.py

    r5755 r5790  
    1 """Script for running tsunami inundation scenario for Dampier, WA, Australia.
    2 
    3 Source data such as elevation and boundary data is assumed to be available in
    4 directories specified by project.py
    5 The output sww file is stored in project.output_run_time_dir
     1"""Script for running a tsunami inundation scenario for carnarvon, WA, Australia.
    62
    73The scenario is defined by a triangular mesh created from project.polygon,
    8 the elevation data and a simulated tsunami generated with URS code.
    9 
    10 Ole Nielsen and Duncan Gray, GA - 2005 and Jane Sexton, Nick Bartzis, GA - 2006
     4the elevation data is compiled into a pts file through build_carnarvon.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)
     8       pts file (build_carnarvon.py)
     9       information from project file
     10Outputs: sww file stored in project.output_run_time_dir
     11The export_results_all.py and get_timeseries.py is reliant
     12on the outputs of this script
     13
     14Ole Nielsen and Duncan Gray, GA - 2005, Jane Sexton, Nick Bartzis, GA - 2006
     15Ole Nielsen, Jane Sexton and Kristy Van Putten - 2008
    1116"""
    1217
    … …  
    3237from Numeric import allclose
    3338from anuga.shallow_water.data_manager import export_grid, create_sts_boundary
    34 
    3539from anuga.pmesh.mesh_interface import create_mesh_from_regions
    3640from anuga.shallow_water.data_manager import start_screen_catcher, copy_code_files,store_parameters
    37 #from anuga_parallel.parallel_api import distribute, numprocs, myid, barrier
    3841from anuga_parallel.parallel_abstraction import get_processor_name
    3942from anuga.caching import myhash
    … …  
    4245from anuga.utilities.polygon import read_polygon, plot_polygons, polygon_area, is_inside_polygon
    4346from anuga.geospatial_data.geospatial_data import find_optimal_smoothing_parameter
    44 
     47from polygon import Polygon_function
     48   
    4549# Application specific imports
    46 import project                 # Definition of file names and polygons
     50import project  # Definition of file names and polygons
    4751numprocs = 1
    4852myid = 0
    … …  
    5054def run_model(**kwargs):
    5155   
    52 
    5356    #------------------------------------------------------------------------------
    5457    # Copy scripts to time stamped output directory and capture screen
    … …  
    5861
    5962    #copy script must be before screen_catcher
    60     #print kwargs
    6163
    6264    print 'output_dir',kwargs['output_dir']
    63     if myid == 0:
    64         copy_code_files(kwargs['output_dir'],__file__,
    65                  dirname(project.__file__)+sep+ project.__name__+'.py' )
    66 
    67         store_parameters(**kwargs)
    68 
    69    # barrier()
     65   
     66    copy_code_files(kwargs['output_dir'],__file__,
     67             dirname(project.__file__)+sep+ project.__name__+'.py' )
     68
     69    store_parameters(**kwargs)
    7070
    7171    start_screen_catcher(kwargs['output_dir'], myid, numprocs)
    7272
    7373    print "Processor Name:",get_processor_name()
    74 
    7574   
    7675    #-----------------------------------------------------------------------
    … …  
    7978
    8079    # Read in boundary from ordered sts file
    81     urs_bounding_polygon=create_sts_boundary(os.path.join(project.boundaries_dir,project.scenario_name))
     80    urs_bounding_polygon=create_sts_boundary(os.path.join(project.boundaries_dir_event,project.scenario_name))
    8281
    8382    # Reading the landward defined points, this incorporates the original clipping
    8483    # polygon minus the 100m contour
    85     landward_bounding_polygon = read_polygon(project.boundaries_dir+'landward_bounding_polygon.txt')
     84    landward_bounding_polygon = read_polygon(project.landward_dir)
    8685
    8786    # Combine sts polyline with landward points
    … …  
    9089    # counting segments
    9190    N = len(urs_bounding_polygon)-1
    92     boundary_tags={'back': [N+1,N+2,N+3], 'side': [N,N+4],'ocean': range(N)}
    93 
    94    
     91
     92    # boundary tags refer to project.landward 4 points equals 5 segments start at N
     93     boundary_tags={'back': [N+1,N+2,N+3], 'side': [N,N+4],'ocean': range(N)}
     94
    9595    #--------------------------------------------------------------------------
    96     # Create the triangular mesh based on overall clipping polygon with a
    97     # tagged
     96    # Create the triangular mesh based on overall clipping polygon with a tagged
    9897    # boundary and interior regions defined in project.py along with
    9998    # resolutions (maximal area of per triangle) for each polygon
    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
    104     if myid == 0:
    105    
    106         print 'start create mesh from regions'
    107 
    108         create_mesh_from_regions(bounding_polygon,
    109                              boundary_tags=boundary_tags,
    110                              maximum_triangle_area=project.res_poly_all,
    111                              interior_regions=project.interior_regions,
    112                              filename=project.meshes_dir_name+'.msh',
    113                              use_cache=True,
    114                              verbose=True)
    115    # barrier()
    116 
    117 ##        covariance_value,alpha = find_optimal_smoothing_parameter (data_file= kwargs['elevation_file'],
    118 ##                                alpha_list=[0.001, 0.01, 0.1, 0.15, 0.2, 0.25, 0.3, 0.4, 0.5],
    119 ##                                mesh_file = project.meshes_dir_name+'.msh')
    120 ##        print 'optimal alpha', covariance_value,alpha       
    121 
     103       
     104    print 'start create mesh from regions'
     105
     106    create_mesh_from_regions(bounding_polygon,
     107                         boundary_tags=boundary_tags,
     108                         maximum_triangle_area=project.res_poly_all,
     109                         interior_regions=project.interior_regions,
     110                         filename=project.meshes_dir_name,
     111                         use_cache=True,
     112                         verbose=True)
     113   
    122114    #-------------------------------------------------------------------------
    123115    # Setup computational domain
    … …  
    125117    print 'Setup computational domain'
    126118
    127     domain = Domain(project.meshes_dir_name+'.msh', use_cache=False, verbose=True)
     119    domain = Domain(project.meshes_dir_name, use_cache=False, verbose=True)
    128120    print 'memory usage before del domain',mem_usage()
    129121       
    … …  
    137129    # Setup initial conditions
    138130    #-------------------------------------------------------------------------
    139     if myid == 0:
    140 
    141         print 'Setup initial conditions'
    142 
    143         from polygon import Polygon_function
    144         #following sets the stage/water to be offcoast only
    145         IC = Polygon_function( [(project.poly_mainland, 0)], default = kwargs['tide'],
    146                                  geo_reference = domain.geo_reference)
    147         domain.set_quantity('stage', IC)
    148         #domain.set_quantity('stage',kwargs['tide'] )
    149         domain.set_quantity('friction', kwargs['friction'])
    150        
    151         print 'Start Set quantity',kwargs['elevation_file']
    152 
    153         domain.set_quantity('elevation',
    154                             filename = kwargs['elevation_file'],
    155                             use_cache = False,
    156                             verbose = True,
    157                             alpha = kwargs['alpha'])
    158         print 'Finished Set quantity'
    159     #barrier()
    160 
     131    print 'Setup initial conditions'
     132
     133    # sets the initial stage in the offcoast region only
     134    IC = Polygon_function( [(project.poly_mainland, 0),(project.poly_island1, 0)
     135                            ,(project.poly_island2, 0)], default = kwargs['tide'],
     136                             geo_reference = domain.geo_reference)
     137    domain.set_quantity('stage', IC)
     138    #domain.set_quantity('stage',kwargs['tide'] )
     139    domain.set_quantity('friction', kwargs['friction'])
     140   
     141    print 'Start Set quantity',kwargs['elevation_file']
     142
     143    domain.set_quantity('elevation',
     144                        filename = kwargs['elevation_file'],
     145                        use_cache = False,
     146                        verbose = True,
     147                        alpha = kwargs['alpha'])
     148    print 'Finished Set quantity'
     149
     150##   #------------------------------------------------------
     151##    # Distribute domain to implement parallelism !!!
    161152##    #------------------------------------------------------
    162 ##    # Create x,y,z file of mesh vertex!!!
    163 ##    #------------------------------------------------------
    164 ##        coord = domain.get_vertex_coordinates()
    165 ##        depth = domain.get_quantity('elevation')
    166 ##       
    167 ##        # Write vertex coordinates to file
    168 ##        filename=project.vertex_filename
    169 ##        fid=open(filename,'w')
    170 ##        fid.write('x (m), y (m), z(m)\n')
    171 ##        for i in range(len(coord)):
    172 ##            pt=coord[i]
    173 ##            x=pt[0]
    174 ##            y=pt[1]
    175 ##            z=depth[i]
    176 ##            fid.write('%.6f,%.6f,%.6f\n' %(x, y, z))
    177153##
    178 
    179     #------------------------------------------------------
    180     # Distribute domain to implement parallelism !!!
    181     #------------------------------------------------------
    182 
    183     if numprocs > 1:
    184         domain=distribute(domain)
     154##    if numprocs > 1:
     155##        domain=distribute(domain)
    185156
    186157    #------------------------------------------------------
    … …  
    188159    #------------------------------------------------------
    189160    print 'domain id', id(domain)
    190     domain.set_name(kwargs['aa_scenario_name'])
     161    domain.set_name(kwargs['scenario_name'])
    191162    domain.set_datadir(kwargs['output_dir'])
    192     domain.set_default_order(2) # Apply second order scheme
    193     domain.set_minimum_storable_height(0.01) # Don't store anything less than 1cm
     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
    194165    domain.set_store_vertices_uniquely(False)
    195166    domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
    … …  
    203174    print 'domain id', id(domain)
    204175   
    205     boundary_urs_out=project.boundaries_dir_name
     176    boundary_urs_out=project.boundaries_dir_event + sep + project.scenario_name
    206177
    207178    Br = Reflective_boundary(domain)
    … …  
    235206        domain.write_boundary_statistics(tags = 'ocean')
    236207
     208    # these outputs should be checked with the resultant inundation map
    237209    x, y = domain.get_maximum_inundation_location()
    238210    q = domain.get_maximum_inundation_elevation()
    239 
    240211    print 'Maximum runup observed at (%.2f, %.2f) with elevation %.2f' %(x,y,q)
    241212
    242     print 'That took %.2f seconds' %(time.time()-t0)
     213    print 'Simulation took %.2f seconds' %(time.time()-t0)
    243214
    244215    #kwargs 'completed' must be added to write the final parameters to file
    245216    kwargs['completed']=str(time.time()-t0)
    246    
    247     if myid==0:
    248         store_parameters(**kwargs)
    249    # barrier
    250    
     217     
     218    store_parameters(**kwargs)
     219     
    251220    print 'memory usage before del domain1',mem_usage()
    252221   
    … …  
    256225   
    257226    kwargs={}
    258     kwargs['est_num_trigs']=project.trigs_min
    259     kwargs['num_cpu']=numprocs
    260     kwargs['host']=project.host
    261     kwargs['res_factor']=project.res_factor
    262     kwargs['starttime']=project.starttime
    263     kwargs['yieldstep']=project.yieldstep
    264227    kwargs['finaltime']=project.finaltime
    265    
    266228    kwargs['output_dir']=project.output_run_time_dir
    267     kwargs['elevation_file']=project.combined_dir_name+'.txt'
    268     kwargs['boundary_file']=project.boundaries_in_dir_name + '.sww'
    269     kwargs['file_name']=project.home+'detail.csv'
    270     kwargs['aa_scenario_name']=project.scenario_name
    271     kwargs['ab_time']=project.time
    272     kwargs['res_factor']= project.res_factor
     229    kwargs['elevation_file']=project.combined_dir_name+'.pts'
     230    kwargs['scenario_name']=project.scenario_name
    273231    kwargs['tide']=project.tide
    274     kwargs['user']=project.user
    275232    kwargs['alpha'] = project.alpha
    276233    kwargs['friction']=project.friction
    277     kwargs['time_thinning'] = project.time_thinning
    278     kwargs['dir_comment']=project.dir_comment
    279     kwargs['export_cellsize']=project.export_cellsize
    280    
    281 
     234    #kwargs['num_cpu']=numprocs
     235    #kwargs['host']=project.host
     236    #kwargs['starttime']=project.starttime
     237    #kwargs['yieldstep']=project.yieldstep
     238    #kwargs['user']=project.user
     239    #kwargs['time_thinning'] = project.time_thinning
     240     
    282241    run_model(**kwargs)
    283242     
    284     #barrier
     243   
Note: See TracChangeset for help on using the changeset viewer.