Changeset 9332


Ignore:
Timestamp:
Sep 11, 2014, 11:44:09 AM (11 years ago)
Author:
davies
Message:

Template scripts now read in special bounding_polygon shapefile format which also gives tags

Location:
trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami1
Files:
6 added
4 edited

Legend:

Unmodified
Added
Removed

  • trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami1/scripts/project.py

    r9331 r9332  
    1414from anuga_parallel import myid, barrier, send, receive, numprocs
    1515# We need the boundary tags info to correctly cache the mesh
    16 from setup_boundary_conditions import boundary_info
     16#from setup_boundary_conditions import boundary_info
     17import user_functions
    1718
    1819
    … …  
    3435## Total model runtime in seconds (starts at time=0.)
    3536#
    36 finaltime=0.01*3600.0
     37finaltime=0.05*3600.0
    3738
    3839## Model time between output file timesteps in seconds
    … …  
    5455## Bounding polygon with default mesh resolution (as a max triangle area in m^2)
    5556#
    56 # Filename can be .shp or .csv
    57 #
    58 bounding_polygon_file='../anuga/polygons/bounding_polygon.csv'
     57# Filename MUST be a .shp file, and contain a set of lines making up the boundary polygon
     58# The lines each have an attribute giving the boundary tag, in the attribute table
     59# column named (the value of boundary_tags_attribute_name)
     60#
     61# The starting/ending points on each line MUST EXACTLY AGREE so the whole
     62# thing makes a polygon.
     63#
     64# Make it in GIS by first making the bounding polygon, then splitting it at
     65# the appropriate locations, then adding the boundary tag attribute.
     66#
     67#bounding_polygon_file='../anuga/polygons/bounding_polygon.csv'
     68bounding_polygon_and_tags_file='../anuga/boundaries/Patong_bounding_polygon/Patong_bounding_polygon.shp'
     69boundary_tags_attribute_name='bndryTag'
     70
    5971default_res=0.5*500**2 
    6072
    … …  
    188200
    189201
    190 bounding_polygon=su.read_polygon(bounding_polygon_file)
     202#bounding_polygon=su.read_polygon(bounding_polygon_file)
     203bounding_polygon_and_tags =\
     204    user_functions.read_boundary_tags_line_shapefile(bounding_polygon_and_tags_file, boundary_tags_attribute_name)
     205bounding_polygon = bounding_polygon_and_tags[0]
     206boundary_tags = bounding_polygon_and_tags[1]
     207
    191208breaklines=su.readListOfBreakLines(breakline_files)
    192209riverwalls, riverwall_par=su.readListOfRiverWalls(riverwall_csv_files) #{} # Dictionary
    … …  
    218235mesh_id_hash=hashlib.md5(json.dumps(\
    219236            [bounding_polygon, interior_regions, riverwalls, \
    220              breaklines, region_point_areas, default_res, boundary_info]) ).hexdigest()
     237             breaklines, region_point_areas, default_res, boundary_tags]) ).hexdigest()
    221238
    222239# Fix the output tif bounding polygon

  • trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami1/scripts/setup_boundary_conditions.py

    r9329 r9332  
    22
    33 THE USER MUST EDIT THIS FILE TO SET THE BOUNDARY CONDITIONS
    4 
    5  Both the boundary_info list and the function 'setup_boundary_conditions'
    6  require adjustment
    7 
    8  Specify boundary conditions
    94
    105"""
    … …  
    138import anuga
    149import anuga.utilities.spatialInputUtil as su
    15 
    16 #
    17 # For each boundary type, we have a list with a
    18 #  1) Boundary tag
    19 #  3) Points -- for each point, the closest segment on the bounding_polygon
    20 #      will be assigned to the boundary type. The default boundary type has
    21 #     'Default' in place of points
    22 #
    23 
    24 #@@@@@@@@@@@@@@@@@#
    25 #START EDITING HERE
    26 #@@@@@@@@@@@@@@@@@#
    27 
    28 ## OCEAN BOUNDARY
    29 # These points define the 'ocean_boundary': 1 point per edge on the bounding polygon
    30 #
    31 ocean_boundary_points_file='../anuga/boundaries/ocean_boundary_points/ocean_boundary_points.shp'
    32 ocean_boundary_points=su.read_points(ocean_boundary_points_file)
    33 
    34 
    35 
    36 ## BOUNDARY TAGS
    37 # Polygon edges will be determined from the points
    38 # Any boundary polygon edge which is not the closest edge to a point will be
    39 # assigned to the 'Default' boundary
    40 #
    41 boundary_info=[
    42               ['ocean_waterlevel_boundary', ocean_boundary_points],
    43               ['vertices_without_tagpoints', 'Default']
    44              ]
    4510
    4611
    … …  
    5318
    5419        2) Next make sure the 'boundary_tags_and_conditions' dictionary
    55            associates the correct 'tags' (defined in boundary_info) to the
    56            correct boundary condition type
     20           associates the correct 'tags' (contained as attribute names in the
     21           bounding_polygon_and_tags_file) with the correct boundary condition type
    5722    """
     23
     24#@@@@@@@@@@@@@@@@@#
     25#START EDITING HERE
     26#@@@@@@@@@@@@@@@@@#
    5827
    5928    # Example:
    6029    # "Standard" boundary conditions
    6130    Br = anuga.Reflective_boundary(domain)
    62     Bs = anuga.Transmissive_stage_zero_momentum_boundary(domain) # Be careful with this one
     31    #Bs = anuga.Transmissive_stage_zero_momentum_boundary(domain) # Be careful with this one
    6332   
    6433    # Example:
    6534    # Transmissive with stage being set by a function 'stage_function'
    6635    def stage_function(t):
    67         return min( max( t-120., 0.)/360., 1.)
     36        return min( max( t-60., 0.)/60., 1.)
    6837    B_t_stage_function = anuga.shallow_water.boundaries.Transmissive_n_momentum_zero_t_momentum_set_stage_boundary(domain, stage_function)
    6938
    7039    ## Here we associate each tag in boundary_info with its boundary condition defined above
    7140    boundary_tags_and_conditions={
    72                                   'ocean_waterlevel_boundary': B_t_stage_function,
    73                                   'vertices_without_tagpoints':Br
     41                                  'Ocean': B_t_stage_function,
     42                                  'Br':Br
    7443                                  }
    7544
    … …  
    7948
    8049    ## Check that all tags in boundary_info have been set
    81     for i in range(len(boundary_info)):
    82         tag=boundary_info[i][0]
     50    for i in range(len(project.boundary_tags)):
     51        tag=project.boundary_tags.keys()[i]
    8352        if not boundary_tags_and_conditions.has_key(tag):
    8453            msg = 'Need to set boundary_tags_and_conditions for tag = ' + tag

  • trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami1/scripts/setup_mesh.py

    r9329 r9332  
    3939#####################################################################################
    4040
    41 def make_boundary_tags(bounding_polygon, setup_boundary_conditions):
    42     """
    43         Make the boundary tags
    44     """
    45 
    46     # Close the bounding polygon
    47     bounding_polygon_closed=bounding_polygon+[bounding_polygon[0]]
    48 
    49     # Assign all tags to 'Default', but remove them from this list during
    50     # assignment
    51     default_tags=range(len(bounding_polygon))
    52 
    53     # Dict for boundary tags
    54     boundary_tags={}
    55 
    56     # Information on the tags
    57     boundary_info=setup_boundary_conditions.boundary_info
    58 
    59     for i in range(len(boundary_info)):
    60         bi=boundary_info[i]
    61         # Shorthand for i'th boundary tag
    62         tag_i=bi[0]
    63         # Shorthand for i'th set of boundary points
    64         PTS_i=bi[1]
    65 
    66         # Add the boundary tag to the list
    67         boundary_tags[tag_i]=[]
    68 
    69         if(PTS_i=='Default' or i==(len(boundary_info)-1)):
    70             # Assign default boundaries
    71             assert (i==(len(boundary_info)-1) and PTS_i=='Default'),\
    72                 'Default boundary condition must be at the end of the boundary_info list'
    73             boundary_tags[tag_i]=default_tags
    74         else:
    75             # Loop over all points in the boundary condition
    76             for j in range(len(PTS_i)):
    77                 new_index=su.find_nearest_segment(PTS_i[j], bounding_polygon_closed)[1]
    78                 boundary_tags[tag_i]=boundary_tags[tag_i]+[new_index]     
    79                 default_tags.remove(new_index)
    80 
    81     return boundary_tags   
     41#def make_boundary_tags(bounding_polygon, setup_boundary_conditions):
     42#    """
     43#        Make the boundary tags
     44#    """
     45#
     46#    # Close the bounding polygon
     47#    bounding_polygon_closed=bounding_polygon+[bounding_polygon[0]]
     48#
     49#    # Assign all tags to 'Default', but remove them from this list during
     50#    # assignment
     51#    default_tags=range(len(bounding_polygon))
     52#
     53#    # Dict for boundary tags
     54#    boundary_tags={}
     55#
     56#    # Information on the tags
     57#    boundary_info=setup_boundary_conditions.boundary_info
     58#
     59#    for i in range(len(boundary_info)):
     60#        bi=boundary_info[i]
     61#        # Shorthand for i'th boundary tag
     62#        tag_i=bi[0]
     63#        # Shorthand for i'th set of boundary points
     64#        PTS_i=bi[1]
     65#
     66#        # Add the boundary tag to the list
     67#        boundary_tags[tag_i]=[]
     68#
     69#        if(PTS_i=='Default' or i==(len(boundary_info)-1)):
     70#            # Assign default boundaries
     71#            assert (i==(len(boundary_info)-1) and PTS_i=='Default'),\
     72#                'Default boundary condition must be at the end of the boundary_info list'
     73#            boundary_tags[tag_i]=default_tags
     74#        else:
     75#            # Loop over all points in the boundary condition
     76#            for j in range(len(PTS_i)):
     77#                new_index=su.find_nearest_segment(PTS_i[j], bounding_polygon_closed)[1]
     78#                boundary_tags[tag_i]=boundary_tags[tag_i]+[new_index]     
     79#                default_tags.remove(new_index)
     80#
     81#    return boundary_tags   
    8282
    8383#################################################################################################
    … …  
    9191
    9292    # Create the boundary tags
    93     boundary_tags=make_boundary_tags(project.bounding_polygon, setup_boundary_conditions)
     93    #boundary_tags=make_boundary_tags(project.bounding_polygon, setup_boundary_conditions)
    9494   
    9595    # Ensure mesh_breaklines include riverwalls and breaklines
    … …  
    9898    ## Make the mesh
    9999    anuga.create_mesh_from_regions(project.bounding_polygon,
    100                              boundary_tags=boundary_tags,
     100                             boundary_tags=project.boundary_tags,
    101101                             maximum_triangle_area=project.default_res,
    102102                             filename=project.meshname,

  • trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami1/scripts/user_functions.py

    r9329 r9332  
    3939    # Make a newline
    4040    if(myid==0): print ''
     41
     42############################################################
     43
     44def read_boundary_tags_line_shapefile(shapefile_name, tag_attribute='bndryTag'):
     45    """
     46    Read in the boundary lines + tags from an appropriately structured line shapefile
     47
     48    Return the bounding polygon + dict with tags + indices
     49    """
     50
     51    ## Step 1: Read the data
     52    import ogr
     53   
     54    data_source=ogr.Open(shapefile_name)
     55
     56    assert data_source.GetLayerCount()==1, ' Bounding polygon + boundary tags file can only have 1 layer'
     57
     58    layer=data_source.GetLayer()
     59
     60    boundary_line_and_tags=[]
     61
     62    for feature in layer:
     63        line = feature.GetGeometryRef().GetPoints()
     64        line = [ list(l) for l in line]
     65        tag = feature.GetField(tag_attribute)
     66        boundary_line_and_tags.append([line, tag])
     67
     68
     69    ## Step 2: Convert to bounding polygon + boundary tags in ANUGA format
     70    boundary_segnum = len(boundary_line_and_tags)
     71    # Initial values
     72    bounding_polygon=boundary_line_and_tags[0][0]
     73    boundary_tags={boundary_line_and_tags[0][1]: range(len(bounding_polygon)-1) }
     74    #
     75    #import pdb
     76    #pdb.set_trace()
     77    for i in range(boundary_segnum-1):
     78        # Find the line which starts at the end point of bounding_polygon
     79        for j in range(boundary_segnum):
     80            blj = boundary_line_and_tags[j][0]
     81            if blj[0]==bounding_polygon[-1]:
     82                # Append all but the first point of blj to bounding_polygon
     83                old_lbp=len(bounding_polygon)
     84                extra_polygon = [blj[k] for k in range(1, len(blj))]
     85                bounding_polygon = bounding_polygon + extra_polygon
     86                # If we are on the last segment, then drop the final point
     87                # (for consistency with ANUGA format)
     88                if i==(boundary_segnum-2):
     89                    bounding_polygon = [bounding_polygon[k] for k in range(len(bounding_polygon)-1)]
     90                    new_tags = range(old_lbp-1, len(bounding_polygon))
     91                else:
     92                    new_tags = range(old_lbp-1, len(bounding_polygon)-1)
     93               
     94                # Either the key already exists and we append edge indices, OR
     95                # we add the key and the edge indices
     96                key = boundary_line_and_tags[j][1]
     97                if boundary_tags.has_key(key):
     98                    boundary_tags[key] = boundary_tags[key]+ new_tags
     99                else:
     100                    boundary_tags[key] = new_tags
     101
     102    return bounding_polygon, boundary_tags
Note: See TracChangeset for help on using the changeset viewer.