Changeset 3788


Ignore:
Timestamp:
Oct 16, 2006, 2:12:04 PM (18 years ago)
Author:
sexton
Message:

path fixups

Location:
anuga_work/production
Files:
7 edited

Legend:

Unmodified
Added
Removed

  • anuga_work/production/broome_2006/project.py

    r3769 r3788  
     1# -*- coding: cp1252 -*-
    12"""Common filenames and locations for topographic data, meshes and outputs.
    23"""
    34
    45from os import sep, environ, getenv, getcwd
    5 from os.path import expanduser, basename
    6 #from anuga.utilities.polygon import read_polygon
     6from os.path import expanduser
    77import sys
    8 from anuga.pmesh.create_mesh import convert_from_latlon_to_utm
    9 from anuga.coordinate_transforms.redfearn import degminsec2decimal_degrees
    10 from time import localtime, strftime
    11 from anuga.geospatial_data.geospatial_data import *
     8from time import localtime, strftime, gmtime
     9from anuga.utilities.polygon import read_polygon, plot_polygons, polygon_area, is_inside_polygon
     10from anuga.coordinate_transforms.redfearn import degminsec2decimal_degrees, convert_points_from_latlon_to_utm
     11
     12if sys.platform == 'win32':
     13    home = getenv('INUNDATIONHOME')
     14    user = getenv('USERPROFILE')
     15
     16else:   
     17    home = getenv('INUNDATIONHOME', sep+'d'+sep+'xrd'+sep+'gem'+sep+'2'+sep+'ramp'+sep+'risk_assessment_methods_project'+sep+'inundation')     
     18    user = getenv('LOGNAME')
     19    print 'USER:', user
     20
     21# INUNDATIONHOME is the inundation directory, not the data directory.
     22home += sep +'data'
    1223
    1324#Making assumptions about the location of scenario data
    … …  
    1526scenario_dir_name = 'broome_tsunami_scenario_2006'
    1627
    17 # all data to be delivered by National Mapping
    18 # onshore data from 30m DTED level 2
    19 onshore_name_dted = 'broome_onshore_30m_dted' 
    20 onshore_name_dli = 'broome_onshore_20m_dli'
     28# onshore data provided by WA DLI
     29onshore_name = '' # original
    2130
    22 # offshore data from GA digitised charts
    23 offshore_name1 = 'broome_offshore_points'
     31# offshore data provided by WA DPI
     32offshore_name_dpi1 = ''
    2433
    25 # offshore data from AHO fairsheets
    26 offshore_name2 = 'broome_offshore_points_fairsheet'
     34# AHO data
     35offshore_name1 = 'xy100003760'
    2736
    28 # coastline developed from aerial photography and 1.5m DLI contour
    29 coast_name = 'broome_coastline_points'
     37# developed by NM&I
     38coast_name = 'coastline_points'
    3039
    31 boundary_basename = 'SU-AU_clip'
     40boundary_basename = 'SU-AU' # Mw ?
    3241
    3342#swollen/ all data output
    … …  
    3544codename = 'project.py'
    3645
    37 if sys.platform == 'win32':
    38     home = getenv('INUNDATIONHOME') #Sandpit's parent dir     
    39     user = getenv('USERPROFILE')
    40 else:
    41     # update to perlite 2
    42     home = getenv('INUNDATIONHOME', sep+'d'+sep+'cit'+sep+'2'+sep+'cit'+sep+'inundation'+sep+'data')     
    43     user = getenv('LOGNAME')
    44 
    45 # INUNDATIONHOME is the inundation directory, not the data directory.
    46 home += sep +'data'
    47 
    4846#Derive subdirectories and filenames
    49 time = strftime('%Y%m%d_%H%M%S',localtime()) #gets time for new dir
    50 outputtimedir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'outputs'+sep+time+sep
    51 
     47local_time = strftime('%Y%m%d_%H%M%S',gmtime())
    5248meshdir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'meshes'+sep
    5349datadir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'topographies'+sep
    5450gaugedir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'gauges'+sep
    5551polygondir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'polygons'+sep
    56 boundarydir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'boundaries'+sep
     52boundarydir = home+sep+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'boundaries'+sep
    5753outputdir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'outputs'+sep
    58 tidedir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'tide_data'+sep
     54outputtimedir = outputdir + local_time + sep
     55polygondir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'polygons'+sep
    5956
    60 gauge_filename = gaugedir + 'gauge_location_broome.csv'
    61 buildings_filename = gaugedir + 'broome_res.csv'
    62 buildings_filename_out = 'broome_res_modified.csv'
    63 buildings_filename_damage_out = 'broome_res_modified_damage.csv'
    64 community_filename = gaugedir + 'CHINS_v2.csv'
    65 community_scenario = gaugedir + 'community_pt_hedland.csv'
    66 tidal_filename = tidedir + 'pt_hedland_tide.txt'
     57gauge_filename = gaugedir + 'broome_gauges.csv'
    6758
    68 meshname = meshdir + basename
    69 onshore_dem_name = datadir + onshore_name_dli
    70 offshore_dem_name1 = datadir + offshore_name1
    71 offshore_dem_name2 = datadir + offshore_name2
     59codedir = getcwd()+sep                           
     60codedirname = codedir + 'project.py'
     61meshname = outputtimedir + 'mesh_' + basename
     62
     63# Necessary if using point datasets, rather than grid
     64onshore_dem_name = datadir + onshore_name
     65offshore_dem_name_local1 = datadir + offshore_name_dpi1
     66offshore_dem_name_aho1 = datadir + offshore_name1
     67
    7268coast_dem_name = datadir + coast_name
    73 combined_dem_name = datadir + 'broome_combined_elevation'
    74 outputname = outputtimedir + basename  #Used by post processing
    7569
    76 # for ferret2sww
    77 south = degminsec2decimal_degrees(-20,30,0)
    78 north = degminsec2decimal_degrees(-17,10,0)
    79 west = degminsec2decimal_degrees(117,00,0)
    80 east = degminsec2decimal_degrees(120,00,0)
     70combined_dem_name   = datadir + 'broome_combined_elevation'
    8171
    82 # region to export (used from export_results.py)
    83 e_min_area = 648000
    84 e_max_area = 675000
    85 n_min_area = 7745000
    86 n_max_area = 7761000
     72###############################
     73# Domain definitions
     74###############################
    8775
    88 export_region = [[e_min_area,n_min_area],[e_min_area,n_max_area],[e_max_area,n_max_area],[e_max_area,n_min_area]]
    89                  
    90 refzone = 50
     76# bounding box for clipping MOST/URS output (much bigger than study area)
     77south = degminsec2decimal_degrees(-19,0,0)
     78north = degminsec2decimal_degrees(-17,15,0)
     79west  = degminsec2decimal_degrees(120,0,0)
     80east  = degminsec2decimal_degrees(122,30,0)
    9181
    92 from anuga.coordinate_transforms.redfearn import redfearn
    93 # boundary up to 50 m contour
    94 lat1_50 = degminsec2decimal_degrees(-19,20,0)
    95 lat2_50 = degminsec2decimal_degrees(-19,30,0)
    96 lat3_50 = degminsec2decimal_degrees(-19,45,0)
    97 lon1_50 = degminsec2decimal_degrees(119,05,0)
    98 lon2_50 = degminsec2decimal_degrees(118,20,0)
    99 lon3_50 = degminsec2decimal_degrees(117,45,0)
    100 z, easting, northing = redfearn(lat1_50, lon1_50)
    101 d0_50 = [easting, northing]
    102 z, easting, northing = redfearn(lat2_50, lon2_50)
    103 d1_50 = [easting, northing]
    104 z, easting, northing= redfearn(lat3_50, lon3_50)
    105 d2_50 = [easting, northing]
     82d0 = [south, west]
     83d1 = [south, east]
     84d2 = [north, east]
     85d3 = [north, west]
     86poly_bc, zone = convert_points_from_latlon_to_utm([d0, d1, d2, d3])
     87refzone = zone
    10688
    107 d4_50 = [285000, 7585000]
    108 d6_50 = [330000, 7605000]
    109 #bounding_poly50 = [p0_50, p1_50, p2_50, d6_50, d5, d4_50]
     89# bounding polygon for study area
     90polyAll = read_polygon(polygondir+'extent.csv')
    11091
    111 d0 = [763852.0, 7934358.0]
    112 d1 = [710987.0, 7925797.0]
    113 d2 = [658264.0, 7926314.0]
    114 d3 = [552686.0, 7871580.0]
    115 #d4 = [604415.81, 7733013.56]
    116 d4 = [638000.0, 7733013.56]
    117 #d5 = [656561.15, 7732615.11]
    118 d5 = [662000.0, 7732615.11]
    119 #d6 = [708940.32, 7750510.33]
    120 d6 = [690000.0, 7740510.33]
    121 #polyAll = [d0, d1, d2, d3, d4, d5, d6]
    122 #polyAll = [d0_50, d1_50, d2_50, d4, d5, d6]
    123 # from Hamish
    124 h0=[629262.17, 7747205.47]
    125 h1=[552686.00, 7871579.99] #d3
    126 h2=[658264.00, 7926314.00] #d2
    127 h3=[710986.99, 7925796.99] #d1
    128 h4=[763851.99, 7934357.99] #d0
    129 h5=[701485.21, 7770656.86]
    130 h6=[698273.75, 7762227.38]
    131 h7=[698194.23, 7762018.65]
    132 h8=[691627.41, 7744781.98]
    133 h9=[679220.75, 7743604.59]
    134 h10=[653512.59, 7740528.56]
    135 h11=[634777.71, 7738247.17]
    136 h12=[629443.86, 7746910.37]
    137 h13=[629396.84, 7746986.75]
    138 h14=[629352.32, 7747059.06]
    139 h15=[629276.24, 7747182.63]
    140 h16=[629262.17, 7747205.47] #repeat of h0
    141 # using Hamish's new bounding polygon
    142 #polyAll = [d0_50, d1_50, d2_50, h16,h15,h14,h13,h12,h11,h10,h9,h8,h7,h6,h5]
    143 polyAll = [d0_50, d1_50, d2_50, h16,h11,h8,h6, h5]
     92# plot bounding polygon and make sure BC info surrounds it
     93plot_polygons([polyAll, poly_bc],'boundingpoly',verbose=False)
     94print 'Area of bounding polygon', polygon_area(polyAll)/1000000.0
    14495
    145 #Interior region - Broome centre
    146 i0 = [668000, 7757000]
    147 i1 = [659000, 7755000]
    148 i2 = [660000, 7749000]
    149 i3 = [667000, 7746000]
    150 i4 = [678000, 7751000]
     96###################################################################
     97# Clipping regions for export to asc and regions for clipping data
     98###################################################################
    15199
    152 poly_broome = [i0, i1, i2, i3, i4]
     100# clipping 20m onshore data set
     101#eastingmin = 520000
     102#eastingmax = 536000
     103#northingmin = 5245000
     104#northingmax = 5260000
    153105
    154 #Are there other significant features?
    155 j0 = [670000, 7760000]
    156 j1 = [633000, 7745000]
    157 j2 = [665000, 7743000]
    158 j3 = [690000, 7755000]
     106###############################
     107# Interior region definitions
     108###############################
    159109
    160 poly_region = [j0, j1, j2, j3]
     110# broome digitized polygons
     111poly_broome1 = read_polygon(polygondir+'Broome_Local_Polygon.csv')
     112poly_broome2 = read_polygon(polygondir+'Broome_Close2.csv')
     113poly_broome3 = read_polygon(polygondir+'Broome_Coast.csv')
     114poly_broome4 = read_polygon(polygondir+'Cable_Beach.csv')
     115
     116plot_polygons([polyAll,poly_broome1,poly_broome2,poly_broome3,poly_broome4],'boundingpoly2',verbose=False)
     117print 'Area of local polygon', polygon_area(poly_broome1)/1000000.0
     118print 'Area of close polygon', polygon_area(poly_broome2)/1000000.0
     119print 'Area of coastal polygon', polygon_area(poly_broome3)/1000000.0
     120print 'Area of cable beach polygon', polygon_area(poly_broome4)/1000000.0
     121
     122for i in poly_broome3:
     123    v = is_inside_polygon(i,poly_broome1, verbose=False)
     124    if v == False: print v

  • anuga_work/production/broome_2006/run_broome.py

    r3535 r3788  
    66
    77The scenario is defined by a triangular mesh created from project.polygon,
    8 the elevation data and a simulated submarine landslide.
     8the elevation data and a tsunami wave generated by MOST.
    99
    1010Ole Nielsen and Duncan Gray, GA - 2005 and Nick Bartzis, GA - 2006
    1111"""
    12 #-------------------------------------------------------------------------------# Import necessary modules
     12def number_mesh_triangles(interior_regions, bounding_poly, remainder_res):
     13    from anuga.utilities.polygon import polygon_area
     14   
     15    # check if any of the regions fall inside one another
     16    no_triangles = 0.0
     17    area = polygon_area(bounding_poly)
     18    for i,j in interior_regions:
     19        this_area = polygon_area(i)
     20        no_triangles += this_area/j
     21        area -= this_area
     22        print j, this_area/1000000., area/1000000.
     23    no_triangles += area/remainder_res
     24    return int(no_triangles/0.7)
     25#-------------------------------------------------------------------------------
     26# Import necessary modules
    1327#-------------------------------------------------------------------------------
    1428
    1529# Standard modules
    16 from os import sep
    17 from os.path import dirname, basename
     30import os
    1831import time
    19 
    20 # Related major packages
    21 from anuga.pyvolution.shallow_water import Domain, Reflective_boundary, \
    22                             Dirichlet_boundary, Time_boundary, File_boundary
    23 from anuga.pyvolution.data_manager import convert_dem_from_ascii2netcdf, dem2pts
    24 from anuga.pyvolution.combine_pts import combine_rectangular_points_files
    25 from anuga.pyvolution.pmesh2domain import pmesh_to_domain_instance
    2632from shutil import copy
    2733from os import mkdir, access, F_OK
     34import sys
     35
     36# Related major packages
     37from anuga.shallow_water import Domain, Reflective_boundary, \
     38                            Dirichlet_boundary, Time_boundary, File_boundary
     39from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf, dem2pts
     40from anuga.abstract_2d_finite_volumes.combine_pts import combine_rectangular_points_files
    2841from anuga.geospatial_data.geospatial_data import *
    29 import sys
    30 from anuga.pyvolution.util import Screen_Catcher
     42from anuga.abstract_2d_finite_volumes.util import Screen_Catcher
    3143
    3244# Application specific imports
    … …  
    4153if access(project.outputtimedir,F_OK) == 0 :
    4254    mkdir (project.outputtimedir)
    43 copy (dirname(project.__file__) +sep+ project.__name__+'.py', project.outputtimedir + project.__name__+'.py')
    44 copy (__file__, project.outputtimedir + basename(__file__))
    45 print 'project.outputtimedir',project.outputtimedir
    46 
    47 # normal screen output is stored in
     55copy (project.codedirname, project.outputtimedir + project.codename)
     56copy (project.codedir + 'run_broome.py', project.outputtimedir + 'run_broome.py')
     57print'output dir', project.outputtimedir
     58
     59#normal screen output is stored in
    4860screen_output_name = project.outputtimedir + "screen_output.txt"
    4961screen_error_name = project.outputtimedir + "screen_error.txt"
    5062
    51 # used to catch screen output to file
     63#used to catch screen output to file
    5264sys.stdout = Screen_Catcher(screen_output_name)
    5365sys.stderr = Screen_Catcher(screen_error_name)
     66
    5467print 'USER:    ', project.user
    5568
    … …  
    5871#
    5972# Convert ASC 2 DEM 2 PTS using source data and store result in source data
    60 # Do for coarse and fine data
    61 # Fine pts file to be clipped to area of interest
    6273#-------------------------------------------------------------------------------
    6374
    6475# filenames
     76#onshore_dem_name = project.onshore_dem_name
     77#coast_points = project.coast_dem_name
     78#offshore_points = project.offshore_dem_name
    6579meshname = project.meshname+'.msh'
    66 source_dir = project.boundarydir
    67 
    68 # creates DEM from asc data
    69 convert_dem_from_ascii2netcdf(project.onshore_dem_name, use_cache=True, verbose=True)
    70 
    71 #creates pts file from DEM
    72 dem2pts(project.onshore_dem_name,
     80#source_dir = project.boundarydir
     81
     82copied_files = False
     83'''
     84# creates DEM from asc data
     85convert_dem_from_ascii2netcdf(onshore_dem_name, use_cache=True, verbose=True)
     86
     87#creates pts file for onshore DEM
     88dem2pts(onshore_dem_name,
    7389        easting_min=project.eastingmin,
    7490        easting_max=project.eastingmax,
    7591        northing_min=project.northingmin,
    7692        northing_max= project.northingmax,
    77         use_cache=True,
     93        use_cache=True, 
    7894        verbose=True)
    7995
    80 print 'create G1'
    81 G1 = Geospatial_data(file_name = project.offshore_dem_name1 + '.xya')
    82 print 'create G2'
    83 G2 = Geospatial_data(file_name = project.offshore_dem_name2 + '.xya')
    84 print 'create G3'
    85 G3 = Geospatial_data(file_name = project.onshore_dem_name + '.pts')
    86 print 'create G4'
    87 G4 = Geospatial_data(file_name = project.coast_dem_name + '.xya')
    88 print 'add G1+G2+G3+G4'
     96print 'create offshore'
     97G1 = Geospatial_data(file_name = project.offshore_dem_name + '.xya')
     98print 'create onshore'
     99G2 = Geospatial_data(file_name = project.onshore_dem_name + '.pts')
     100print 'create coast'
     101G3 = Geospatial_data(file_name = project.coast_dem_name + '.xya')
     102print 'add'
    89103G = G1 + G2 + G3 + G4
    90 print 'export G'
     104print 'export points'
    91105G.export_points_file(project.combined_dem_name + '.pts')
    92 
    93 #-------------------------------------------------------------------------------                                 
     106'''
     107#----------------------------------------------------------------------------
    94108# Create the triangular mesh based on overall clipping polygon with a tagged
    95109# boundary and interior regions defined in project.py along with
    … …  
    98112
    99113from anuga.pmesh.mesh_interface import create_mesh_from_regions
    100 
    101 region_res = 500000
     114remainder_res = 100000
     115local_res = 15000
     116broome_res = 2500
    102117coast_res = 500
    103 broome_res = 5000
    104 interior_regions = [[project.poly_broome, broome_res],
    105                     [project.poly_region, region_res]]
     118beach_res = 500
     119interior_regions = [[project.poly_broome1, local_res],
     120                    [project.poly_broome2, broome_res],
     121                    [project.poly_broome3, coast_res],
     122                    [project.poly_broome4, beach_res]]
    106123
    107124print 'number of interior regions', len(interior_regions)
    108 
    109 from anuga.utilities.polygon import plot_polygons
    110 if sys.platform == 'win32':
    111     #figname = project.outputtimedir + 'pt_hedland_polys'
    112     figname = 'broome_polys_test'
    113     plot_polygons([project.polyAll,project.poly_broome,project.poly_region],
    114               figname,
    115               verbose = True)   
    116 
    117 print 'start create mesh from regions'
     125print 'estimate of number of triangles', \
     126      number_mesh_triangles(interior_regions, project.polyAll, remainder_res)
     127
     128import sys; sys.exit()
     129
    118130from caching import cache
    119131_ = cache(create_mesh_from_regions,
    120132          project.polyAll,
    121           {'boundary_tags': {'topright': [0], 'topleft': [1],
    122                              'left': [2], 'bottom0': [3],
    123                              'bottom1': [4], 'bottom2': [5],
    124                              'bottom3': [6], 'right': [7]},
    125            'maximum_triangle_area': 250000,
    126            'filename': meshname,           
     133           {'boundary_tags': {'e0': [0], 'e1': [1], 'e2': [2],
     134                              'e3': [3], 'e4':[4], 'e5': [5],
     135                              'e6': [6]},
     136           'maximum_triangle_area': remainder_res,
     137           'filename': meshname,
    127138           'interior_regions': interior_regions},
    128           verbose = True, evaluate=True)
    129 
     139          verbose = True, evaluate=False)
     140print 'created mesh'
     141import sys; sys.exit()
    130142#-------------------------------------------------------------------------------                                 
    131143# Setup computational domain
    132144#-------------------------------------------------------------------------------                                 
    133 domain = Domain(meshname, use_cache = False, verbose = True)
    134 
    135 print domain.statistics()
     145domain = Domain(meshname, use_cache = True, verbose = True)
     146
    136147print 'Number of triangles = ', len(domain)
    137148print 'The extent is ', domain.get_extent()
    … …  
    141152domain.set_datadir(project.outputtimedir)
    142153domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
     154domain.set_minimum_storable_height(0.01)
    143155
    144156#-------------------------------------------------------------------------------                                 
    … …  
    146158#-------------------------------------------------------------------------------
    147159
    148 tide = 0.
    149 #high
    150 #tide =
    151 #low
    152 #tide =
    153 
     160tide = 0.0
    154161domain.set_quantity('stage', tide)
    155162domain.set_quantity('friction', 0.0)
    156 print 'hi and file',project.combined_dem_name + '.pts'
    157 
    158 domain.set_quantity('elevation',
    159                     filename = project.combined_dem_name + '.pts',
     163domain.set_quantity('elevation', 0.0,
     164                    #filename = project.combined_dem_name + '.pts',
    160165                    use_cache = True,
    161166                    verbose = True,
    … …  
    164169
    165170#-------------------------------------------------------------------------------                                 
    166 # Setup boundary conditions 
    167 #-------------------------------------------------------------------------------
    168 print 'start ferret2sww'
    169 # skipped as results in file SU-AU_clipped is correct for all WA
    170 
    171 from anuga.pyvolution.data_manager import ferret2sww
     171# Setup boundary conditions
     172#-------------------------------------------------------------------------------
     173'''
     174print 'start urs2sww'
     175print '', project.boundary_basename
     176from anuga.shallow_water.data_manager import urs2sww
    172177
    173178south = project.south
    174179north = project.north
    175 west = project.west
    176 east = project.east
     180west  = project.west
     181east  = project.east
    177182
    178183#note only need to do when an SWW file for the MOST boundary doesn't exist
    179 cache(ferret2sww,
     184cache(urs2sww,
    180185      (source_dir + project.boundary_basename,
    181        source_dir + project.boundary_basename+'_'+project.basename),
     186       source_dir + project.boundary_basename),
    182187      {'verbose': True,
    183188       'minlat': south,
    … …  
    185190       'minlon': west,
    186191       'maxlon': east,
    187 #       'origin': project.mesh_origin,
    188        'origin': domain.geo_reference.get_origin(),
     192       #'origin': domain.geo_reference.get_origin(),
    189193       'mean_stage': tide,
    190194       'zscale': 1,                 #Enhance tsunami
    191195       'fail_on_NaN': False,
    192196       'inverted_bathymetry': True},
    193        evaluate = True,
     197      #evaluate = True,
    194198       verbose = True,
    195       dependencies = source_dir + project.boundary_basename + '.sww')
    196 
     199       dependencies = source_dir + project.boundary_basename + '.sww')
     200
     201'''
    197202print 'Available boundary tags', domain.get_boundary_tags()
    198203
    199 Bf = File_boundary(source_dir + project.boundary_basename + '.sww',
    200                     domain, verbose = True)
     204#Bf = File_boundary(source_dir + project.boundary_basename + '.sww',
     205#                    domain, verbose = True)
    201206Br = Reflective_boundary(domain)
    202207Bd = Dirichlet_boundary([tide,0,0])
    203 domain.set_boundary( {'topright': Bf,'topleft': Bf, 'left':  Bd, 'bottom0': Bd,
    204                       'bottom1': Bd, 'bottom2': Bd, 'bottom3': Bd,
    205                         'right': Bd})
     208
     209# 7 min square wave starting at 1 min, 6m high
     210Bw = Time_boundary(domain = domain,
     211                   f=lambda t: [(60<t<480)*10, 0, 0])
     212
     213domain.set_boundary( {'e0': Bd,  'e1': Bd, 'e2': Bd, 'e3': Bd, 'e4': Bd,
     214                      'e5': Bd,  'e6': Bd} )
     215
    206216
    207217#-------------------------------------------------------------------------------                                 
    … …  
    211221t0 = time.time()
    212222
    213 for t in domain.evolve(yieldstep = 240, finaltime = 10800):
     223for t in domain.evolve(yieldstep = 240, finaltime = 6800):
    214224    domain.write_time()
    215     domain.write_boundary_statistics(tags = 'topright')     
    216 
    217 for t in domain.evolve(yieldstep = 120, finaltime = 16200
     225    domain.write_boundary_statistics(tags = 'e14')     
     226
     227for t in domain.evolve(yieldstep = 30, finaltime = 9000
    218228                       ,skip_initial_step = True):
    219229    domain.write_time()
    220     domain.write_boundary_statistics(tags = 'topright')     
    221 
    222 for t in domain.evolve(yieldstep = 60, finaltime = 21600
     230    domain.write_boundary_statistics(tags = 'e14')     
     231
     232for t in domain.evolve(yieldstep = 240, finaltime = 15000
    223233                       ,skip_initial_step = True):
    224234    domain.write_time()
    225     domain.write_boundary_statistics(tags = 'topright')   
     235    domain.write_boundary_statistics(tags = 'e14')
    226236   
    227 for t in domain.evolve(yieldstep = 120, finaltime = 27000
    228                        ,skip_initial_step = True):
    229     domain.write_time()
    230     domain.write_boundary_statistics(tags = 'topright')     
    231 
    232 for t in domain.evolve(yieldstep = 240, finaltime = 36000
    233                        ,skip_initial_step = True):
    234     domain.write_time()
    235     domain.write_boundary_statistics(tags = 'topright')   
    236  
    237237print 'That took %.2f seconds' %(time.time()-t0)
    238238

  • anuga_work/production/onslow_2006/project.py

    r3769 r3788  
    55from os import sep, environ, getenv, getcwd
    66from os.path import expanduser
    7 #from anuga.utilities.polygon import read_polygon
     7from anuga.utilities.polygon import read_polygon, polygon_area
    88import sys
    99
    … …  
    163163
    164164polyAll = [d0, d1, d2, d3, d4, d5, d6]
    165 
     165print 'bounding polygon area', polygon_area(polyAll)/1000000.0
    166166polygons = [polyAll, export_region]
    167167figname = 'checking.png'
    … …  
    208208
    209209poly_onslow = [i0, i1, i2, i3, i4, i5]
    210 
     210print 'onslow polygon area', polygon_area(poly_onslow)/1000000.0
    211211#Thevenard Island
    212212j0 = [294000, 7629000]
    … …  
    216216
    217217poly_thevenard = [j0, j1, j2, j3]
    218 
     218print 'thevenard polygon area', polygon_area(poly_thevenard)/1000000.0
    219219'''
    220220# Direction Is
    … …  
    237237#poly_coast = [l0, l1, l2, l3, l4, l5]
    238238poly_coast = [l0, l1, l2, l3, l4]
    239 
     239print 'coast polygon area', polygon_area(poly_coast)/1000000.0
    240240#general coast and local area to onslow region
    241241m0 = [270000, 7581000]
    … …  
    247247
    248248poly_region = [m0, m1, m2, m3, m4, m5]
     249print 'region polygon area', polygon_area(poly_region)/1000000.0

  • anuga_work/production/onslow_2006/run_onslow.py

    r3535 r3788  
    2323
    2424# Related major packages
    25 from anuga.pyvolution.shallow_water import Domain, Reflective_boundary, \
     25from anuga.shallow_water import Domain, Reflective_boundary, \
    2626                            Dirichlet_boundary, Time_boundary, File_boundary
    27 from anuga.pyvolution.data_manager import convert_dem_from_ascii2netcdf, dem2pts
    28 from anuga.pyvolution.combine_pts import combine_rectangular_points_files
    29 from anuga.pyvolution.pmesh2domain import pmesh_to_domain_instance
     27from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf, dem2pts
     28from anuga.abstract_2d_finite_volumes.combine_pts import combine_rectangular_points_files
    3029from anuga.geospatial_data.geospatial_data import *
    31 from anuga.pyvolution.util import Screen_Catcher
     30from anuga.abstract_2d_finite_volumes.util import Screen_Catcher
    3231
    3332# Application specific imports
    … …  
    144143#-------------------------------------------------------------------------------                                 
    145144
    146 domain = pmesh_to_domain_instance(meshname, Domain,
    147                                   use_cache = False,
    148                                   verbose = True)
     145#domain = pmesh_to_domain_instance(meshname, Domain,
     146#                                  use_cache = False,
     147#                                  verbose = True)
     148
     149domain = Domain(meshname, use_cache = False, verbose = True)
    149150
    150151print 'Number of triangles = ', len(domain)

  • anuga_work/production/pt_hedland_2006/project.py

    r3770 r3788  
    77#from anuga.utilities.polygon import read_polygon
    88import sys
    9 from anuga.pmesh.create_mesh import convert_from_latlon_to_utm
     9from anuga.coordinate_transforms.redfearn import convert_points_from_latlon_to_utm
    1010from anuga.coordinate_transforms.redfearn import degminsec2decimal_degrees
    1111from time import localtime, strftime

  • anuga_work/production/pt_hedland_2006/run_pt_hedland.py

    r3535 r3788  
    1919
    2020# Related major packages
    21 from anuga.pyvolution.shallow_water import Domain, Reflective_boundary, \
     21from anuga.shallow_water import Domain, Reflective_boundary, \
    2222                            Dirichlet_boundary, Time_boundary, File_boundary
    23 from anuga.pyvolution.data_manager import convert_dem_from_ascii2netcdf, \
     23from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf, \
    2424     dem2pts
    25 from anuga.pyvolution.combine_pts import combine_rectangular_points_files
    26 from anuga.pyvolution.pmesh2domain import pmesh_to_domain_instance
     25from anuga.abstract_2d_finite_volumes.combine_pts import combine_rectangular_points_files
    2726from shutil import copy
    2827from os import mkdir, access, F_OK
    2928from anuga.geospatial_data.geospatial_data import *
    3029import sys
    31 from anuga.pyvolution.util import Screen_Catcher
     30from anuga.abstract_2d_finite_volumes.util import Screen_Catcher
    3231
    3332# Application specific imports

  • anuga_work/production/sydney_2006/project.py

    r3769 r3788  
    88import sys
    99
    10 from anuga.pmesh.create_mesh import convert_from_latlon_to_utm
     10from anuga.coordinate_transforms.redfearn import degminsec2decimal_degrees, convert_from_latlon_to_utm
    1111               
    1212
Note: See TracChangeset for help on using the changeset viewer.