source: inundation/pmesh/create_mesh.py @ 2299

"""Create mesh for lwru2 validation of the Oshiri island tsunami
"""

#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
# Assume that the root AnuGA dir (inundation) is included in your pythonpath
#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

#FIXME (DSG): move the coordinate transforms package to the root AnuGA dir 
from coordinate_transforms.geo_reference import Geo_reference
from coordinate_transforms.redfearn import redfearn


def convert_points_from_latlon_to_utm(polygon, refzone):
    points = []
    for point in polygon:
        zone, easting, northing  = redfearn(point[0], point[1])
        #FIXME: Use point.latitude etc once we have a proper point set
        assert zone == refzone
        points.append([easting, northing])

    return points    





def create_region(polygon, tags, refzone):
    """Create dictionary representation of region bounded by given polygon for use with pmesh
    """


    #
    dict = {}
    dict['points'] = polygon

    #Create segments
    #E.g. [[0,1], [1,2], [2,3], [3,0]]
    segments = []
    N = len(polygon)
    for i in range(N):
        lo = i
        hi = (lo + 1) % N
        segments.append( [lo, hi] ) 
        
        
    dict['segments'] = segments


    #Create tags
    #E.g. ['wall', 'wall', 'ocean', 'wall']

    if tags is not None:
        segment_tags = [0]*N
        for key in tags:
            indices = tags[key]
            for i in indices:
                segment_tags[i] = key
    
        dict['segment_tags'] = segment_tags



    return dict



def create_mesh_from_regions(bounding_polygon,
                             boundary_tags,
                             resolution,
                             filename = None,
                             interior_regions = None):
    """Create mesh from bounding polygon, tags for all segments and resolution.

    Polygon is a list of points in latitude and longitudes - decimal degrees

    Boundary tags is a dictionary of symbolic tags. For every tag there is a list of
    indices referring to segments associated with that tag

    Resolution is the maximal area per triangle for the bounding polygon
    (excluding interior regions, see later)

    Interior_regions is a list of tuples consisting of (polygon, resolution)
    for each region to be separately refined.

    This function does not allow segments to share points - use underlying
    pmesh functionality for that

    FIXME - FUTURE. Something like this.
      Use class Point_set

      Accept deg, min, sec, e.g.
      [ [(-20,30,55), (116, 20, 00)], ...]



    
    """

    #from pmesh.mesh import Mesh
    from mesh import Mesh
    from coordinate_transforms.redfearn import redfearn
    from pyvolution.util import populate_polygon
    

    #Make georef
    #FIXME: Pass in geo or create automatically somehow
    import project
    refzone = project.refzone  #FIXME
    mesh_origin = project.mesh_origin
    geo = Geo_reference(xllcorner = mesh_origin[1],  #From dem
                        yllcorner = mesh_origin[2],
                        zone = refzone)



    
    print "***********************"
    print "geo ref", geo
    print "***********************"
    
    m = Mesh(geo_reference=geo)



    #Convert to UTM
    bounding_polygon = convert_points_from_latlon_to_utm(bounding_polygon, refzone)    
    
    
    #Do bounding polygon
    region_dict = create_region(bounding_polygon, boundary_tags, refzone)    
    m.addVertsSegs(region_dict)


    #Find one point inside region automatically
    if interior_regions is not None:
        excluded_polygons = []        
        for P, res in interior_regions:
            polygon = convert_points_from_latlon_to_utm(P, refzone)            
            excluded_polygons.append( polygon )
    else:
        excluded_polygons = None

    from Numeric import array
    
    [inner_point] = populate_polygon(bounding_polygon, 1, exclude = excluded_polygons)
    print 'I:', inner_point, resolution
    inner = m.addRegionEN(inner_point[0], inner_point[1])
    inner.setMaxArea(resolution)



    #Do interior regions
    if interior_regions is not None:    
        for P, res in interior_regions:
            polygon = convert_points_from_latlon_to_utm(P, refzone)                        
            region_dict = create_region(polygon, None, refzone)

            m.addVertsSegs(region_dict)

            [inner_point] = populate_polygon(polygon, 1)
            print 'I', inner_point, res
            X = m.addRegionEN(inner_point[0], inner_point[1])
            X.setMaxArea(res)
            

    
    m.generateMesh('pzq28.0za1000000a')
    
    m.export_mesh_file(filename)