source: inundation-numpy-branch/pmesh/create_mesh.py @ 3330

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

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

from coordinate_transforms.redfearn import redfearn

def convert_points_from_latlon_to_utm(polygon, refzone=None):
    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

        if refzone is None:
            refzone = zone
        else:
            assert zone == refzone
        
        points.append([easting, northing])

    return points, zone    




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,
                             UTM = False,
                             refzone = 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 mesh import Mesh
    from coordinate_transforms.redfearn import redfearn
    from utilities.polygon import populate_polygon
    from utilities.numerical_tools import ensure_numeric
    from coordinate_transforms.geo_reference import Geo_reference


    bounding_polygon = ensure_numeric(bounding_polygon)
    #print 'refzone', refzone

    #Convert to UTM
    if not UTM:
        bounding_polygon, refzone = convert_points_from_latlon_to_utm(bounding_polygon)
    else:
       assert refzone is not None, 'Refzone must be specified, got %s' %refzone  


    #FIXME: Must give georeference here.
    #Why can't Mesh work out the ll corners?

    xllcorner = min(bounding_polygon[:,0])
    yllcorner = min(bounding_polygon[:,1])    
    
    geo = Geo_reference(xllcorner = xllcorner,
                        yllcorner = yllcorner,
                        zone = refzone)

    print 'geo reference derived from bounding polygon', geo
    m = Mesh(geo_reference=geo)
    
    
    
    #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:
            if not UTM:            
                polygon, _ = convert_points_from_latlon_to_utm(P, refzone)
            else:
                polygon = P
            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:
            if not UTM:                        
                polygon, _ = convert_points_from_latlon_to_utm(P, refzone)
            else:
                polygon = P
                
            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)