Changeset 6279


Ignore:
Timestamp:
Feb 5, 2009, 2:52:42 PM (15 years ago)
Author:
kristy
Message:

Updated standardised, not tested yet

Location:
anuga_work/production/busselton/standardised_version
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • anuga_work/production/busselton/standardised_version/build_busselton.py

    r6276 r6279  
    3030# output to file
    3131#------------------------------------------------------------------------------
    32 copy_code_files(project.output_build_time_dir,__file__,
     32copy_code_files(project.output_build,__file__,
    3333                os.path.dirname(project.__file__)+os.sep+\
    3434                project.__name__+'.py' )
    35 start_screen_catcher(project.output_build_time_dir)
     35start_screen_catcher(project.output_build)
    3636
    3737
     
    4444#------------------------------------------------------------------------------
    4545
    46 print 'project.poly_all', project.poly_all
    47 print 'project.combined_dir_name', project.combined_dir_name
     46print 'project.bounding_polygon', project.bounding_polygon
     47print 'project.combined_elevation_basename', project.combined_elevation_basename
    4848
    4949# Create Geospatial data from ASCII files
    5050geospatial_data = {}
    5151for filename in project.ascii_grid_filenames:
    52     absolute_filename = project.topographies_dir + filename
     52    absolute_filename = project.topographies_folder + os.sep + filename
    5353    convert_dem_from_ascii2netcdf(absolute_filename,
    5454                                  basename_out=absolute_filename,
     
    6262# Create Geospatial data from TXT files
    6363for filename in project.point_filenames:
    64     absolute_filename = project.topographies_dir + filename
     64    absolute_filename = project.topographies_folder  + os.sep + filename
    6565    geospatial_data[filename] = Geospatial_data(file_name=absolute_filename,
    6666                                                verbose=True)
     
    7171#-------------------------------------------------------------------------------
    7272
    73 print 'Add geospatial objects except', project.offshore_name5
    74 G = None
    75 for key in geospatial_data:
    76     if key != project.offshore_name5:
    77         G += geospatial_data[key]
     73print 'Add geospatial objects' # except', project.offshore_name5
     74G = geospatial_data
     75##None
     76##for key in geospatial_data:
     77##    if key != project.offshore_name5:
     78##        G += geospatial_data[key]
    7879       
    7980print 'Clip combined geospatial data'
    80 G_clip = G.clip_outside(project.poly_aoi1)
    81 G_all = G_clip + geospatial_data[project.offshore_name5]   # G_off5
    82 G_clipped = G_all.clip(project.poly_all)
     81##G_clip = G.clip_outside(project.poly_aoi1)
     82##G_all = G_clip + geospatial_data[project.offshore_name5]   # G_off5
     83#G_clipped = G_all.clip(project.poly_all)
     84G_clip = G.clip(project.bounding_polygon)
     85
    8386
    8487print 'Export combined DEM file'
    85 G_clipped.export_points_file(project.combined_dir_name + '.pts')
     88G_clip.export_points_file(project.combined_dir_name + '.pts')
    8689print 'Do txt version too'
    8790# Use for comparision in ARC
    88 G_clipped.export_points_file(project.combined_dir_name + '.txt')
     91G_clip.export_points_file(project.combined_dir_name + '.txt')
    8992
  • anuga_work/production/busselton/standardised_version/project.py

    r6278 r6279  
    2929run_time = time+'_run'
    3030
    31 #------------------------------------------------------------------------------
    32 # Initial Conditions
    33 #------------------------------------------------------------------------------
    34 
    3531# this section needs to be updated to reflect the modelled community.
    3632# Note, the user needs to set up the directory system accordingly
     
    3935scenario_folder = 'busselton_tsunami_scenario'
    4036
     37#------------------------------------------------------------------------------
     38# Initial Conditions
     39#------------------------------------------------------------------------------
    4140# Model specific parameters. One or all can be changed each time the
    4241# run_scenario script is executed
     
    5049finaltime=80000         # final time for simulation
    5150
    52 setup='trial'  # Final can be replaced with trial or basic.
     51setup='final'  # Final can be replaced with trial or basic.
    5352               # Either will result in a coarser mesh that will allow a
    5453               # faster, but less accurate, simulation.
     
    5655if setup =='trial':
    5756    print'trial'
    58     res_factor=100
     57    scale_factor=100
    5958    time_thinning=96
    6059    yieldstep=240
    6160if setup =='basic':
    6261    print'basic'
    63     res_factor=4
     62    scale_factor=4
    6463    time_thinning=12
    6564    yieldstep=120
    6665if setup =='final':
    6766    print'final'
    68     res_factor=1
     67    scale_factor=1
    6968    time_thinning=4
    7069    yieldstep=60
     
    7675# Important to distinguish each run - ensure str(user) is included!
    7776# Note, the user is free to include as many parameters as desired
    78 output_comment= join('_', setup, '_', str(tide), '_', str(event_number),
    79                      '_', str(user)
     77output_comment= ('_' + setup + '_' + str(tide)+ '_' + str(event_number) +
     78                 '_' + str(user))
    8079
    8180#------------------------------------------------------------------------------
    8281# Input Data
    8382#------------------------------------------------------------------------------
    84 
    85 # elevation data - used in build_busselton.py
     83# ELEVATION DATA
     84# Used in build_busselton.py
    8685# Format for ascii grids, as produced in ArcGIS + a projection file
    87 ascii_grid_filenames = [onshore_name,   # Topo
    88                         offshore_name5] # Busselton Topo
     86ascii_grid_filenames = ['busselton_v2',   # Topo
     87                        'grid_250m_clip'] # Busselton Topo
    8988
    9089# Format for point is x,y,elevation (with header)
    91 point_filenames = [coast_name,     # Coastline
    92                    coast_name1,    # Beach survey
    93                    offshore_name,  # Bathymetry
    94                    offshore_name1, # Bathymetry Charts
    95                    offshore_name2, # Digitised Fairsheet
    96                    offshore_name3, # 250m
    97                    offshore_name4] # Bunbury DPI
    98 
    99 # Land used to set the initial stage/water to be offcoast only
     90point_filenames = ['Busselton_Contour0.txt',     # Coastline
     91                   'Busselton_BeachSurvey.txt',    # Beach survey
     92                   'Busselton_NavyFinal.txt',  # Bathymetry
     93                   'Busselton_Chart.txt', # Bathymetry Charts
     94                   'Busselton_Digitised.txt', # Digitised Fairsheet
     95                   'Busselton_250m.txt', # 250m
     96                   'Bunbury_TIN.txt', # Bunbury aoi TIN'd in ArcGIS
     97                   'Busselton_TIN.txt', # Busselton aoi TIN'd in ArcGIS
     98                   'XYAHD_clip.txt'] # To extend boundary
     99
     100
     101# LAND - used to set the initial stage/water to be offcoast only
    100102# Used in run_busselton,py
    101103# Format for points easting,northing (no header)
    102 mainland_polygon_filename = 'initial_condition.csv'
    103 land_polygon_filename = 'initial_condition_land.csv'
    104 
    105 # Bounding polygon for data clipping and estimate of triangles in mesh
     104land_initial_conditions_filename = [['initial_condition_extend.csv', 0],
     105                                    ['initial_condition_marina.csv', 0]]
     106
     107# BOUNDING POLYGON - for data clipping and estimate of triangles in mesh
    106108# Used in build_busselton.py
    107109# Format for points easting,northing (no header)
    108110bounding_polygon_filename = 'bounding_polygon.csv'
    109111
    110 # Interior regions, for designing the mesh - used in run_busselton.py
     112# INTERIOR REGIONS -  for designing the mesh
     113# Used in run_busselton.py
    111114# Format for points easting,northing (no header)                   
    112115interior_regions_data = [['busselton_1km.csv', 500],
     
    116119                         ['island1.csv', 10000],
    117120                         ['island2.csv', 10000],
    118                          ['coast_5km_d20m.csv', 40000
    119 # gauges - used in get_timeseries.py
     121                         ['coast_5km_d20m.csv', 40000]]
     122
     123# GAUGES - for creating timeseries at a specific point
     124# Used in get_timeseries.py
    120125# Format easting,northing,name,elevation (with header)
    121126gauges_filename = 'gauges.csv'
    122127
    123 # buildings - used in run_building_inundation.py
    124 # Format latitude,longitude etc (not geographic)
     128# BUILDINGS EXPOSURE - for identifying inundated houses
     129# Used in run_building_inundation.py
     130# Format latitude,longitude etc (geographic)
    125131building_exposure_filename = 'busselton_res_clip.csv' #from NEXIS
    126132
    127 # BOUNDING POLYGON - used in build_boundary.py and run_busselton.py respectively
     133# BOUNDING POLYGON
     134# used in build_boundary.py and run_busselton.py respectively
    128135# NOTE: when files are put together the points must be in sequence
    129136# For ease go clockwise!
     
    132139# thinned ordering file from Hazard Map (geographic)
    133140# Format is index,latitude,longitude (with header)
    134 urs_order_filename = 'thinned_boundary_ordering.csv'
     141urs_order_filename = 'thinned_boundary_ordering_extend.csv'
    135142
    136143#landward bounding points
    137144# Format easting,northing (no header)
    138 landward_boundary_filename = 'landward_boundary.csv'
     145landward_boundary_filename = 'landward_boundary_extend.csv'
     146
     147#------------------------------------------------------------------------------
     148# Clipping regions for export to asc and regions for clipping data
     149# Final inundation maps should only be created in regions of the finest mesh
     150#------------------------------------------------------------------------------
     151
     152# ASCII export grid for Busselton
     153xminBusselton = 340000
     154xmaxBusselton = 352000
     155yminBusselton = 6271500
     156ymaxBusselton = 6280000
     157
     158# ASCII export grid for Bunbury
     159xminBunbury = 369000
     160xmaxBunbury = 381000
     161yminBunbury = 6308000
     162ymaxBunbury = 6316500
    139163
    140164
     
    149173# Directory Structure
    150174#------------------------------------------------------------------------------
    151 anuga_folder = join(home, state, scenario, 'anuga')
     175anuga_folder = join(home, state, scenario_folder, 'anuga')
    152176topographies_folder = join(anuga_folder, 'topographies')
    153177polygons_folder = join(anuga_folder, 'polygons')
     
    179203# The absolute pathname for the output folder names
    180204# Used for build_busselton.py
    181 output_build = join(output_folder, build_time, dir_comment)
     205output_build = output_folder + sep + build_time + output_comment
    182206# Used for run_busselton.py
    183 output_run = join(output_folder, run_time, dir_comment)
     207output_run = output_folder + sep + run_time + output_comment
    184208# Used by post processing
    185209output_run_time = join(output_run, scenario_name)
     
    193217
    194218#------------------------------------------------------------------------------
    195 # Interior region
    196 #------------------------------------------------------------------------------
    197 
    198 # Land, to set the initial stage/water to be offcoast only
    199 mainland_polygon = read_polygon(join(polygons_folder,
    200                                      mainland_polygon_filename))
    201 
    202 # Land, to set the initial stage/water to be offcoast only
    203 land_polygon = read_polygon(join(polygons_folder,
    204                                    land_polygon_filename))
     219# Reading polygons and creating interior regions
     220#------------------------------------------------------------------------------
     221
     222land_initial_conditions = []
     223for filename, MSL in land_initial_conditions_filename:
     224    polygon = read_polygon(join(polygons_folder, filename))
     225    land_initial_conditions.append([filename, MSL])
    205226
    206227# Initial bounding polygon for data clipping
     
    210231
    211232interior_regions = []
    212 for (filename, maxarea) in interior_regions_data:
    213     polygon = read_polygon(join(polygon_folder, filename)
    214     maxarea = maxarea*scale_factor
    215     interior_regions.append([filename, maxarea])
     233for filename, maxarea in interior_regions_data:
     234    polygon = read_polygon(join(polygons_folder, filename))
     235    interior_regions.append([polygon, maxarea*scale_factor])
    216236                     
    217    
    218237trigs_min = number_mesh_triangles(interior_regions,
    219238                                  bounding_polygon, bounding_maxarea)
    220239print 'min estimated number of triangles', trigs_min
    221240   
    222 #------------------------------------------------------------------------------
    223 # Clipping regions for export to asc and regions for clipping data
    224 # Final inundation maps should only be created in regions of the finest mesh
    225 #------------------------------------------------------------------------------
    226 
    227 # ASCII export grid for Busselton
    228 xminBusselton = 340000
    229 xmaxBusselton = 352000
    230 yminBusselton = 6271500
    231 ymaxBusselton = 6280000
    232 
    233 # ASCII export grid for Bunbury
    234 xminBunbury = 369000
    235 xmaxBunbury = 381000
    236 yminBunbury = 6308000
    237 ymaxBunbury = 6316500
    238 
     241
  • anuga_work/production/busselton/standardised_version/run_busselton.py

    r6267 r6279  
    2626# Related major packages
    2727from anuga.interface import create_domain_from_regions
     28from anuga.interface import Transmissive_stage_zero_momentum_boundary
    2829from anuga.interface import Dirichlet_boundary
    2930from anuga.interface import Reflective_boundary
     
    4041
    4142
    42 #-----------------------------------------------------------------------
     43#------------------------------------------------------------------------------
    4344# Copy scripts to time stamped output directory and capture screen
    4445# output to file. Copy script must be before screen_catcher
    45 #-----------------------------------------------------------------------
    46 copy_code_files(project.output_run_time_dir, __file__,
     46#------------------------------------------------------------------------------
     47copy_code_files(project.output_run, __file__,
    4748                os.path.dirname(project.__file__)+os.sep+\
    4849                project.__name__+'.py' )
    49 start_screen_catcher(project.output_run_time_dir, 0, 1)
     50start_screen_catcher(project.output_run, 0, 1)
    5051
    5152
    52 #--------------------------------------------------------------------------
     53#------------------------------------------------------------------------------
    5354# Create the computational domain based on overall clipping polygon with
    5455# a tagged boundary and interior regions defined in project.py along with
    5556# resolutions (maximal area of per triangle) for each polygon
    56 #--------------------------------------------------------------------------
     57#------------------------------------------------------------------------------
    5758print 'Create computational domain'
    5859
    5960# Read in boundary from ordered sts file
    60 urs_bounding_polygon = create_sts_boundary(project.urs_boundary_name)
     61event_sts = create_sts_boundary(project.event_sts)
    6162
    6263# Reading the landward defined points, this incorporates the original clipping
    6364# polygon minus the 100m contour
    64 landward_bounding_polygon = read_polygon(project.landward_dir)
     65landward_boundary = read_polygon(project.landward_boundary)
    6566
    6667# Combine sts polyline with landward points
    67 bounding_polygon = urs_bounding_polygon + landward_bounding_polygon
     68bounding_polygon_sts = event_sts + landward_boundary
    6869
    6970# Number of boundary segments
    70 N = len(urs_bounding_polygon)-1
     71N = len(event_sts)-1
    7172
    7273# Boundary tags refer to project.landward 4 points equals 5 segments start at N
     
    7677
    7778# Build mesh and domain
    78 domain = create_domain_from_regions(bounding_polygon,
     79domain = create_domain_from_regions(bounding_polygon_sts,
    7980                                    boundary_tags=boundary_tags,
    8081                                    maximum_triangle_area=project.res_poly_all,
    8182                                    interior_regions=project.interior_regions,
    82                                     mesh_filename=project.meshes_dir_name,
     83                                    mesh_filename=project.meshes,
    8384                                    use_cache=True,
    8485                                    verbose=True)
     
    8687
    8788domain.set_name(project.scenario_name)
    88 domain.set_datadir(project.output_run_time_dir)
     89domain.set_datadir(project.output_run)
    8990domain.set_minimum_storable_height(0.01)    # Don't store depth less than 1cm
    9091
    9192
    92 #-------------------------------------------------------------------------
     93#------------------------------------------------------------------------------
    9394# Setup initial conditions
    94 #-------------------------------------------------------------------------
     95#------------------------------------------------------------------------------
    9596print 'Setup initial conditions'
    9697
    9798# Set the initial stage in the offcoast region only
    98 IC = Polygon_function([(project.poly_mainland, 0),
    99                        (project.poly_marina, 0)],
     99IC = Polygon_function(project.land_initial_conditions,
    100100                      default=project.tide,
    101101                      geo_reference=domain.geo_reference)
     
    103103domain.set_quantity('friction', project.friction)
    104104domain.set_quantity('elevation',
    105                     filename=project.combined_dir_name+'.pts',
     105                    filename=project.combined_elevation+'.pts',
    106106                    use_cache=True,
    107107                    verbose=True,
     
    109109
    110110
    111 #-------------------------------------------------------------------------
     111#------------------------------------------------------------------------------
    112112# Setup boundary conditions
    113 #-------------------------------------------------------------------------
     113#------------------------------------------------------------------------------
    114114print 'Set boundary - available tags:', domain.get_boundary_tags()
    115115
    116116Br = Reflective_boundary(domain)
    117 Bd = Dirichlet_boundary([project.tide,0,0])
    118 Bf = Field_boundary(project.urs_boundary_name+'.sts',
     117Bt = Transmissive_stage_zero_momentum_boundary(domain)
     118Bf = Field_boundary(project.event_sts+'.sts',
    119119                    domain, mean_stage=project.tide,
    120120                    time_thinning=1,
    121121                    default_boundary=Bd,
    122                     boundary_polygon=bounding_polygon,                   
     122                    boundary_polygon=bounding_polygon_sts,                   
    123123                    use_cache=True,
    124124                    verbose=True)
    125125
    126126domain.set_boundary({'back': Br,
    127                      'side': Bd,
     127                     'side': Bt,
    128128                     'ocean': Bf})
    129129
    130130
    131 #-------------------------------------------------------------------------
     131#------------------------------------------------------------------------------
    132132# Evolve system through time
    133 #-------------------------------------------------------------------------
     133#------------------------------------------------------------------------------
    134134t0 = time.time()
    135135for t in domain.evolve(yieldstep=project.yieldstep,
Note: See TracChangeset for help on using the changeset viewer.