Changeset 2615


Ignore:
Timestamp:
Mar 28, 2006, 11:00:20 AM (18 years ago)
Author:
nick
Message:

update onslow

Location:
production/onslow_2006
Files:
2 added
1 edited
2 moved

Legend:

Unmodified
Added
Removed
  • production/onslow_2006/project_new.py

    r2470 r2615  
    99
    1010from pmesh.create_mesh import convert_points_from_latlon_to_utm
     11
     12from coordinate_transforms.redfearn import degminsec2decimal_degrees
     13
     14from time import localtime, strftime
     15
     16from os import getcwd
    1117               
    1218#Making assumptions about the location of scenario data
    1319scenario_dir_name = 'onslow_tsunami_scenario_2006'
    1420
    15 
    1621# 250m data to be provided
    1722coarsename = 'onsl_bathydem250' # get from Neil/Ingo (DEM or topo data)
    18 """
     23
    1924# 30m data to be provided
    20 finename = 'onshore_30' # get from Neil/Ingo (DEM or topo data)
    21 """
     25onshore_name = 'onslow_onshore_30m_dted' # get from Neil/Ingo (DEM or topo data)
    2226
    23 # clipping region for fine elevation data
    24 west = 240000
    25 east = 340000
    26 south = 7570000
    27 north = 7645000
     27offshore_name = 'onslow_offshore_points'
     28
     29boundary_basename = 'SU-AU'
    2830
    2931#swollen/ all data output
    3032basename = 'source'
     33
     34codename = 'project.py'
    3135
    3236if sys.platform == 'win32':
     
    3640
    3741#Derive subdirectories and filenames
     42timedir = strftime('%Y%m%d_%H%M%S',localtime()) #gets time for new dir
     43outputdir = home+sep+scenario_dir_name+sep+'output'+sep+timedir+sep
    3844meshdir = home+sep+scenario_dir_name+sep+'meshes'+sep
    3945datadir = home+sep+scenario_dir_name+sep+'topographies'+sep
    40 outputdir = home+sep+scenario_dir_name+sep+'output'+sep
     46
    4147polygondir = home+sep+scenario_dir_name+sep+'polygons'+sep
    4248boundarydir = home+sep+scenario_dir_name+sep+'boundaries'+sep
     49
     50codedir = getcwd()+sep
     51                               
     52codedirname = codedir + 'project.py'
    4353
    4454meshname = meshdir + basename
    4555
    4656coarsedemname = datadir + coarsename
    47 """
    48 finedemname = datadir + finename
    49 combineddemname = datadir + 'onslow_combined_elevation'
    50 """
     57
     58onshore_dem_name = datadir + onshore_name
     59
     60offshore_dem_name = datadir + offshore_name
     61
     62combined_dem_name = datadir + 'onslow_combined_elevation'
    5163
    5264outputname = outputdir + basename  #Used by post processing
     
    5466#!gauge_filename = outputdir + 'onslow_gauges.xya'
    5567#!gauge_outname = outputdir + 'gauges_max_output.xya'
     68
     69# clipping region for fine elevation data
     70eastingmin = 250000
     71eastingmax = 330000
     72northingmin = 7580000
     73northingmax = 7635000
     74
     75south = degminsec2decimal_degrees(-22,00,0)
     76north = degminsec2decimal_degrees(-21,10,0)
     77west = degminsec2decimal_degrees(114,30,0)
     78east = degminsec2decimal_degrees(115,30,0)
     79
     80# region for visualisation
     81eminviz = 260000
     82emaxviz = 320000
     83nminviz = 7590000
     84nmaxviz = 7630000
    5685
    5786#Georeferencing
     
    97126poly_direction = [k0, k1, k2, k3]
    98127
    99 #!slump_origin = [385000.0, 6255000.0] #Absolute UTM
  • production/onslow_2006/run_onslow_new.py

    r2470 r2615  
    2121# Related major packages
    2222from pyvolution.shallow_water import Domain, Reflective_boundary, \
    23                             Dirichlet_boundary, Time_boundary
     23                            Dirichlet_boundary, Time_boundary, File_boundary
    2424from pyvolution.data_manager import convert_dem_from_ascii2netcdf, dem2pts
    2525from pyvolution.combine_pts import combine_rectangular_points_files
    2626from pyvolution.pmesh2domain import pmesh_to_domain_instance
     27from geospatial_data import add_points_files
    2728
    2829# Application specific imports
     
    3031from smf import slump_tsunami  # Function for submarine mudslide
    3132
     33from shutil import copy
     34from os import mkdir, access, F_OK
     35
     36from geospatial_data import *
    3237
    3338#-------------------------------------------------------------------------------
     
    4146# filenames
    4247coarsedemname = project.coarsedemname
    43 '''
    44 #finedemname = project.finedemname
    45 '''
     48
     49onshore_dem_name = project.onshore_dem_name
     50
     51offshore_points = project.offshore_dem_name
     52
    4653meshname = project.meshname+'.msh'
    4754
     55source_dir = project.boundarydir
     56
     57# creates copy of code in output dir
     58if access(project.outputdir,F_OK) == 0 :
     59    mkdir (project.outputdir)
     60copy (project.codedirname, project.outputdir + project.codename)
     61copy (project.codedir + 'run_onslow.py', project.outputdir + 'run_onslow.py')
     62
     63
     64'''
    4865# coarse data
    4966convert_dem_from_ascii2netcdf(coarsedemname, use_cache=True, verbose=True)
    5067dem2pts(coarsedemname, use_cache=True, verbose=True)
    51 '''
     68
     69
    5270# fine data (clipping the points file to smaller area)
    53 convert_dem_from_ascii2netcdf(finedemname, use_cache=True, verbose=True)
    54 dem2pts(finedemname,
     71convert_dem_from_ascii2netcdf(onshore_dem_name, use_cache=True, verbose=True)
     72dem2pts(onshore_dem_name,
    5573        easting_min=project.eastingmin,
    5674        easting_max=project.eastingmax,
     
    6078        verbose=True)
    6179
    62 
    63 # combining the coarse and fine data
    64 combine_rectangular_points_files(project.finedemname + '.pts',
    65                                  project.coarsedemname + '.pts',
    66                                  project.combineddemname + '.pts')
    67 '''
    68 
     80        '''
     81print'create G1'
     82G1 = Geospatial_data(file_name = project.offshore_dem_name + '.xya')
     83
     84print'create G2'
     85G2 = Geospatial_data(file_name = project.onshore_dem_name + '.pts')
     86
     87print'add G1+G2'
     88G = G1 + G2
     89
     90print'export G'
     91G.new_export_points_file(project.combined_dem_name + '.pts')
     92
     93'''
     94add_points_files(
     95                  project.offshore_dem_name + '.xya',
     96                  project.onshore_dem_name + '.pts',
     97                  project.combined_dem_name + '.pts')
     98'''
    6999#-------------------------------------------------------------------------------                                 
    70100# Create the triangular mesh based on overall clipping polygon with a tagged
     
    73103#-------------------------------------------------------------------------------
    74104
    75 from pmesh.create_mesh import create_mesh_from_regions
     105from pmesh.mesh_interface import create_mesh_from_regions
    76106
    77107# original
     
    80110                    [project.poly_thevenard, interior_res],
    81111                    [project.poly_direction, interior_res]]
    82 
    83 #FIXME: Fix caching of this one once the mesh_interface is ready
     112                    #[project.testpoly, interior_res]]
     113print 'number of interior regions', len(interior_regions)
     114
    84115from caching import cache
    85116_ = cache(create_mesh_from_regions,
     
    88119                             'left': [2], 'bottom': [3],
    89120                             'bottomright': [4], 'topright': [5]},
    90            'resolution': 100000,
     121           'maximum_triangle_area': 100000,
    91122           'filename': meshname,           
    92            'interior_regions': interior_regions,
    93            'UTM': True,
    94            'refzone': project.refzone},
     123           'interior_regions': interior_regions},
    95124          verbose = True)
    96125
     
    106135print 'Number of triangles = ', len(domain)
    107136print 'The extent is ', domain.get_extent()
     137print domain.statistics()
    108138
    109139domain.set_name(project.basename)
    110140domain.set_datadir(project.outputdir)
    111 domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
     141domain.set_quantities_to_be_stored(['stage'])
    112142
    113143
     
    136166domain.set_quantity('stage', tide)
    137167domain.set_quantity('friction', 0.0)
     168print 'hi and file',project.combined_dem_name + '.pts'
    138169domain.set_quantity('elevation',
    139 #                    filename = project.combineddemname + '.pts',
    140                     filename = project.coarsedemname + '.pts',
    141                     use_cache = True,
    142                     verbose = True
     170#                    0.
     171#                    filename = project.onshore_dem_name + '.pts',
     172                    filename = project.combined_dem_name + '.pts',
     173#                    filename = project.offshore_dem_name + '.pts',
     174                    use_cache = False,
     175                    verbose = True,
     176                    alpha = 0.1
    143177                    )
    144 
     178print 'hi1'
    145179
    146180#-------------------------------------------------------------------------------                                 
     
    156190
    157191cache(ferret2sww,
    158       (source_dir+project.boundary_basename,
    159        project.boundary_basename),
     192      (source_dir + project.boundary_basename,
     193       source_dir + project.boundary_basename),
    160194      {'verbose': True,
    161        'minlat': south-1,
    162        'maxlat': north+1,
    163        'minlon': west-1,
    164        'maxlon': east+1,
    165        'origin': project.mesh_origin,
     195# note didn't work with the below
     196#       'minlat': south - 1,
     197#       'maxlat': north + 1,
     198#       'minlon': west - 1,
     199#       'maxlon': east + 1,
     200       'minlat': south,
     201       'maxlat': north,
     202       'minlon': west,
     203       'maxlon': east,
     204#       'origin': project.mesh_origin,
     205       'origin': domain.geo_reference.get_origin(),
    166206       'mean_stage': tide,
    167207       'zscale': 1,                 #Enhance tsunami
     
    174214print 'Available boundary tags', domain.get_boundary_tags()
    175215
     216Bf = File_boundary(source_dir + project.boundary_basename + '.sww',
     217                    domain, verbose = True)
    176218Br = Reflective_boundary(domain)
    177219Bd = Dirichlet_boundary([tide,0,0])
     
    182224                   f=lambda t: [(60<t<480)*6, 0, 0])
    183225
    184 domain.set_boundary( {'top': Bw, 'topleft': Br,
     226domain.set_boundary( {'top': Bf, 'topleft': Bf,
    185227                             'left': Br, 'bottom': Br,
    186                              'bottomright': Br, 'topright': Br} )
     228                             'bottomright': Br, 'topright': Bf} )
    187229
    188230
     
    190232# Evolve system through time
    191233#-------------------------------------------------------------------------------
    192 
    193234import time
    194235t0 = time.time()
     
    199240   
    200241print 'That took %.2f seconds' %(time.time()-t0)
     242
     243print 'finished'
Note: See TracChangeset for help on using the changeset viewer.