Changeset 3881


Ignore:
Timestamp:
Oct 30, 2006, 9:21:02 AM (18 years ago)
Author:
sexton
Message:

broome updates

Location:
anuga_work/production/broome_2006
Files:
1 added
2 edited

Legend:

Unmodified
Added
Removed
  • anuga_work/production/broome_2006/project.py

    r3788 r3881  
    88from time import localtime, strftime, gmtime
    99from anuga.utilities.polygon import read_polygon, plot_polygons, polygon_area, is_inside_polygon
    10 from anuga.coordinate_transforms.redfearn import degminsec2decimal_degrees, convert_points_from_latlon_to_utm
     10#from anuga.coordinate_transforms.redfearn import degminsec2decimal_degrees, convert_points_from_latlon_to_utm
    1111
    1212if sys.platform == 'win32':
     
    2727
    2828# onshore data provided by WA DLI
    29 onshore_name = '' # original
     29onshore_name = 'dted2_z51' # original
    3030
    31 # offshore data provided by WA DPI
    32 offshore_name_dpi1 = ''
    33 
    34 # AHO data
    35 offshore_name1 = 'xy100003760'
     31# AHO + DPI data
     32offshore_name1 = 'XY100011610'
     33offshore_name2 = 'XY100011611'
     34offshore_name3 = 'XY100011613'
     35offshore_name4 = 'XY100011614'
     36offshore_name5 = 'XY100011616'
     37offshore_name6 = 'XY100011617'
     38offshore_name7 = 'XY100011618'
     39offshore_name8 = 'XY100011621'
     40offshore_name9 = 'XY100011623'
     41offshore_name10 = 'XY100011745'
     42offshore_name11 = 'XY100011746'
     43offshore_name12 = 'XY100017530'
     44offshore_name13 = 'XY100017532'
     45offshore_name14 = 'XY100017538'
     46offshore_name15 = 'XY100017540'
     47offshore_name16 = 'XYBR66'
     48offshore_name17 = 'XYBR70'
     49offshore_name18 = 'XYBR80'
     50offshore_name19 = 'XYBR88'
     51offshore_name20 = 'XYBR93'
     52offshore_name21 = 'XYBR0110'
     53offshore_name22 = 'XYWADPI'
    3654
    3755# developed by NM&I
    38 coast_name = 'coastline_points'
     56coast_name = 'coastline'
    3957
    4058boundary_basename = 'SU-AU' # Mw ?
     
    6381# Necessary if using point datasets, rather than grid
    6482onshore_dem_name = datadir + onshore_name
    65 offshore_dem_name_local1 = datadir + offshore_name_dpi1
    66 offshore_dem_name_aho1 = datadir + offshore_name1
     83offshore_dem_name1 = datadir + offshore_name1
     84offshore_dem_name2 = datadir + offshore_name2
     85offshore_dem_name3 = datadir + offshore_name3
     86offshore_dem_name4 = datadir + offshore_name4
     87offshore_dem_name5 = datadir + offshore_name5
     88offshore_dem_name6 = datadir + offshore_name6
     89offshore_dem_name7 = datadir + offshore_name7
     90offshore_dem_name8 = datadir + offshore_name8
     91offshore_dem_name9 = datadir + offshore_name9
     92offshore_dem_name10 = datadir + offshore_name10
     93offshore_dem_name11 = datadir + offshore_name11
     94offshore_dem_name12 = datadir + offshore_name12
     95offshore_dem_name13 = datadir + offshore_name13
     96offshore_dem_name14 = datadir + offshore_name14
     97offshore_dem_name15 = datadir + offshore_name15
     98offshore_dem_name16 = datadir + offshore_name16
     99offshore_dem_name17 = datadir + offshore_name17
     100offshore_dem_name18 = datadir + offshore_name18
     101offshore_dem_name19 = datadir + offshore_name19
     102offshore_dem_name20 = datadir + offshore_name20
     103offshore_dem_name21 = datadir + offshore_name21
     104offshore_dem_name22 = datadir + offshore_name22
    67105
    68106coast_dem_name = datadir + coast_name
     
    75113
    76114# bounding box for clipping MOST/URS output (much bigger than study area)
    77 south = degminsec2decimal_degrees(-19,0,0)
    78 north = degminsec2decimal_degrees(-17,15,0)
    79 west  = degminsec2decimal_degrees(120,0,0)
    80 east  = degminsec2decimal_degrees(122,30,0)
    81 
    82 d0 = [south, west]
    83 d1 = [south, east]
    84 d2 = [north, east]
    85 d3 = [north, west]
    86 poly_bc, zone = convert_points_from_latlon_to_utm([d0, d1, d2, d3])
    87 refzone = zone
     115##south = degminsec2decimal_degrees(-19,0,0)
     116##north = degminsec2decimal_degrees(-17,15,0)
     117##west  = degminsec2decimal_degrees(120,0,0)
     118##east  = degminsec2decimal_degrees(122,30,0)
     119##
     120##d0 = [south, west]
     121##d1 = [south, east]
     122##d2 = [north, east]
     123##d3 = [north, west]
     124##poly_bc, zone = convert_points_from_latlon_to_utm([d0, d1, d2, d3])
     125##refzone = zone
    88126
    89127# bounding polygon for study area
     
    91129
    92130# plot bounding polygon and make sure BC info surrounds it
    93 plot_polygons([polyAll, poly_bc],'boundingpoly',verbose=False)
     131#plot_polygons([polyAll, poly_bc],'boundingpoly',verbose=False)
    94132print 'Area of bounding polygon', polygon_area(polyAll)/1000000.0
    95133
     
    98136###################################################################
    99137
    100 # clipping 20m onshore data set
    101 #eastingmin = 520000
    102 #eastingmax = 536000
    103 #northingmin = 5245000
    104 #northingmax = 5260000
     138# exporting asc grid
     139eastingmin = 406215.87
     140eastingmax = 440208.78
     141northingmin = 7983427.73
     142northingmax = 8032834.52
     143
     144       
     145       
     146
    105147
    106148###############################
     
    112154poly_broome2 = read_polygon(polygondir+'Broome_Close2.csv')
    113155poly_broome3 = read_polygon(polygondir+'Broome_Coast.csv')
    114 poly_broome4 = read_polygon(polygondir+'Cable_Beach.csv')
     156poly_broome4 = read_polygon(polygondir+'Cable_Beach_revised.csv')
    115157
    116158plot_polygons([polyAll,poly_broome1,poly_broome2,poly_broome3,poly_broome4],'boundingpoly2',verbose=False)
     
    123165    v = is_inside_polygon(i,poly_broome1, verbose=False)
    124166    if v == False: print v
     167
     168def number_mesh_triangles(interior_regions, bounding_poly, remainder_res):
     169    from anuga.utilities.polygon import polygon_area
     170   
     171    # TO DO check if any of the regions fall inside one another
     172    no_triangles = 0.0
     173    area = polygon_area(bounding_poly)
     174    for i,j in interior_regions:
     175        this_area = polygon_area(i)
     176        no_triangles += this_area/j
     177        area -= this_area
     178        print j, this_area/1000000., area/1000000.
     179    no_triangles += area/remainder_res
     180    return int(no_triangles/0.7)
  • anuga_work/production/broome_2006/run_broome.py

    r3788 r3881  
    1010Ole Nielsen and Duncan Gray, GA - 2005 and Nick Bartzis, GA - 2006
    1111"""
    12 def 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)
     12
    2513#-------------------------------------------------------------------------------
    2614# Import necessary modules
     
    7462
    7563# filenames
    76 #onshore_dem_name = project.onshore_dem_name
    77 #coast_points = project.coast_dem_name
    78 #offshore_points = project.offshore_dem_name
     64onshore_dem_name = project.onshore_dem_name
     65coast_points = project.coast_dem_name
    7966meshname = project.meshname+'.msh'
    80 #source_dir = project.boundarydir
    81 
    82 copied_files = False
    83 '''
     67
    8468# creates DEM from asc data
    8569convert_dem_from_ascii2netcdf(onshore_dem_name, use_cache=True, verbose=True)
    8670
    8771#creates pts file for onshore DEM
    88 dem2pts(onshore_dem_name,
    89         easting_min=project.eastingmin,
    90         easting_max=project.eastingmax,
    91         northing_min=project.northingmin,
    92         northing_max= project.northingmax,
    93         use_cache=True,
    94         verbose=True)
     72dem2pts(onshore_dem_name, use_cache=True, verbose=True)
    9573
    9674print 'create offshore'
    97 G1 = Geospatial_data(file_name = project.offshore_dem_name + '.xya')
     75G1 = Geospatial_data(file_name = project.offshore_dem_name1 + '.xya')+\
     76     Geospatial_data(file_name = project.offshore_dem_name2 + '.xya')+\
     77     Geospatial_data(file_name = project.offshore_dem_name3 + '.xya')+\
     78     Geospatial_data(file_name = project.offshore_dem_name4 + '.xya')+\
     79     Geospatial_data(file_name = project.offshore_dem_name5 + '.xya')+\
     80     Geospatial_data(file_name = project.offshore_dem_name6 + '.xya')+\
     81     Geospatial_data(file_name = project.offshore_dem_name7 + '.xya')+\
     82     Geospatial_data(file_name = project.offshore_dem_name8 + '.xya')+\
     83     Geospatial_data(file_name = project.offshore_dem_name9 + '.xya')+\
     84     Geospatial_data(file_name = project.offshore_dem_name10 + '.xya')+\
     85     Geospatial_data(file_name = project.offshore_dem_name11 + '.xya')+\
     86     Geospatial_data(file_name = project.offshore_dem_name12 + '.xya')+\
     87     Geospatial_data(file_name = project.offshore_dem_name13 + '.xya')+\
     88     Geospatial_data(file_name = project.offshore_dem_name14 + '.xya')+\
     89     Geospatial_data(file_name = project.offshore_dem_name15 + '.xya')+\
     90     Geospatial_data(file_name = project.offshore_dem_name16 + '.xya')+\
     91     Geospatial_data(file_name = project.offshore_dem_name17 + '.xya')+\
     92     Geospatial_data(file_name = project.offshore_dem_name18 + '.xya')+\
     93     Geospatial_data(file_name = project.offshore_dem_name19 + '.xya')+\
     94     Geospatial_data(file_name = project.offshore_dem_name20 + '.xya')+\
     95     Geospatial_data(file_name = project.offshore_dem_name21 + '.xya')+\
     96     Geospatial_data(file_name = project.offshore_dem_name22 + '.xya')
    9897print 'create onshore'
    9998G2 = Geospatial_data(file_name = project.onshore_dem_name + '.pts')
     
    101100G3 = Geospatial_data(file_name = project.coast_dem_name + '.xya')
    102101print 'add'
    103 G = G1 + G2 + G3 + G4
     102G = G1 + G2 + G3
    104103print 'export points'
    105104G.export_points_file(project.combined_dem_name + '.pts')
    106 '''
     105
    107106#----------------------------------------------------------------------------
    108107# Create the triangular mesh based on overall clipping polygon with a tagged
     
    112111
    113112from anuga.pmesh.mesh_interface import create_mesh_from_regions
    114 remainder_res = 100000
    115 local_res = 15000
    116 broome_res = 2500
     113remainder_res = 750000
     114local_res = 25000
     115broome_res = 5000
    117116coast_res = 500
    118117beach_res = 500
     
    122121                    [project.poly_broome4, beach_res]]
    123122
    124 print 'number of interior regions', len(interior_regions)
     123from project import number_mesh_triangles
    125124print 'estimate of number of triangles', \
    126125      number_mesh_triangles(interior_regions, project.polyAll, remainder_res)
    127 
    128 import sys; sys.exit()
    129126
    130127from caching import cache
     
    139136          verbose = True, evaluate=False)
    140137print 'created mesh'
    141 import sys; sys.exit()
     138
    142139#-------------------------------------------------------------------------------                                 
    143140# Setup computational domain
     
    161158domain.set_quantity('stage', tide)
    162159domain.set_quantity('friction', 0.0)
    163 domain.set_quantity('elevation', 0.0,
    164                     #filename = project.combined_dem_name + '.pts',
     160domain.set_quantity('elevation',
     161                    filename = project.combined_dem_name + '.pts',
    165162                    use_cache = True,
    166163                    verbose = True,
     
    221218t0 = time.time()
    222219
    223 for t in domain.evolve(yieldstep = 240, finaltime = 6800):
     220for t in domain.evolve(yieldstep = 240, finaltime = 480):
    224221    domain.write_time()
    225222    domain.write_boundary_statistics(tags = 'e14')     
    226 
    227 for t in domain.evolve(yieldstep = 30, finaltime = 9000
    228                        ,skip_initial_step = True):
    229     domain.write_time()
    230     domain.write_boundary_statistics(tags = 'e14')     
    231 
    232 for t in domain.evolve(yieldstep = 240, finaltime = 15000
    233                        ,skip_initial_step = True):
    234     domain.write_time()
    235     domain.write_boundary_statistics(tags = 'e14')
    236223   
    237224print 'That took %.2f seconds' %(time.time()-t0)
Note: See TracChangeset for help on using the changeset viewer.