zoukankan      html  css  js  c++  java
  • 生成凹包

    import time
    from struct import *
    import subprocess
    import fiona
    import math
    import numpy as np
    import pylab as pl
    from osgeo import ogr, gdal
    import shapely.geometry as geometry
    from shapely.geometry import polygon, multipolygon
    from geojson import Feature, Point, FeatureCollection, dumps
    from descartes import PolygonPatch
    from shapely.ops import cascaded_union, polygonize
    from scipy.spatial import Delaunay
    
    #
    input_shapefile = 'E:/test/pointNoRepeat.shp'
    shapefile = fiona.open(input_shapefile)
    points = [geometry.shape(point['geometry']) for point in shapefile]
    x = [p.coords.xy[0] for p in points]
    y = [p.coords.xy[1] for p in points]
    
    
    '''
    point_collection = geometry.MultiPoint(list(points))
    point_collection.envelope
    
    def plot_polygon(polygon):
        fig = pl.figure(figsize=(10,10))
        ax = fig.add_subplot(111)
        margin = .3
        x_min, y_min, x_max, y_max = polygon.bounds
        ax.set_xlim([x_min-margin, x_max+margin])
        ax.set_ylim([y_min-margin, y_max+margin])
        #print(str(x_min) + "  " + str(y_min) + "  " + str(x_max) + "  " + str(y_max))
        #print(polygon)
        patch = PolygonPatch(polygon, fc='#999999', ec='#000000', fill=True, zorder=-1)
        ax.add_patch(patch)
        return fig
    
    plot_polygon(point_collection.envelope)
    pl.plot(x,y,'o', color='#f16824')
    '''
    
    
    #
    def alpha_shape(points, alpha):
        """
        Compute the alpha shape (concave hull) of a set of points.
        @param points: Iterable container of points.
        @param alpha: alpha value to influence the gooeyness of the border. Smaller numbers don't fall inward as much as larger numbers.
            Too large, and you lose everything!
        """
        if len(points) < 4:
            # When you have a triangle, there is no sense in computing an alpha shape.
            return geometry.MultiPoint(list(points)).convex_hull
        def add_edge(edges, edge_points, coords, i, j):
            """
            Add a line between the i-th and j-th points, if not in the list already
            """
            if (i, j) in edges or (j, i) in edges:
                # already added
                return
            edges.add( (i, j) )
            edge_points.append(coords[ [i, j] ])
        coords = np.array([point.coords[0] for point in points])
        #print(coords)
        tri = Delaunay(coords)
        edges = set()
        edge_points = []
        # loop over triangles:
        # ia, ib, ic = indices of corner points of the triangle
        for ia, ib, ic in tri.vertices:
            pa = coords[ia]
            pb = coords[ib]
            pc = coords[ic]
            # Lengths of sides of triangle
            a = math.sqrt((pa[0]-pb[0])**2 + (pa[1]-pb[1])**2)
            b = math.sqrt((pb[0]-pc[0])**2 + (pb[1]-pc[1])**2)
            c = math.sqrt((pc[0]-pa[0])**2 + (pc[1]-pa[1])**2)
            # Semiperimeter of triangle
            s = (a + b + c)/2.0
            # Area of triangle by Heron's formula
            area = math.sqrt(s*(s-a)*(s-b)*(s-c))
            if area <= 0:
                continue;
            circum_r = a*b*c/(4.0*area)
            # Here's the radius filter.
            #print circum_r
            if circum_r < 1.0/alpha:
                add_edge(edges, edge_points, coords, ia, ib)
                add_edge(edges, edge_points, coords, ib, ic)
                add_edge(edges, edge_points, coords, ic, ia)
        m = geometry.MultiLineString(edge_points)
        triangles = list(polygonize(m))
        #
        return cascaded_union(triangles), edge_points
    
    #
    concave_hull, edge_points = alpha_shape(points, alpha=0.38)
    #
    featCollection = FeatureCollection([])
    featCollection["crs"] = {"type":"name","properties":{"name":"EPSG:2385"}}
    if hasattr(concave_hull, "geoms"):
        polys = concave_hull.geoms
        fig = pl.figure(figsize=(10,10))
        ax = fig.add_subplot(111)
        margin = .3
        if len(concave_hull.bounds)>0:
            #
            x_min, y_min, x_max, y_max = concave_hull.bounds
            ax.set_xlim([x_min-margin, x_max+margin])
            ax.set_ylim([y_min-margin, y_max+margin])
            for poly in polys:
                patch = PolygonPatch(poly, fc='#999999', ec='#000000', fill=True, zorder=-1)
                ax.add_patch(patch)
                #
                geojs_geom = poly.__geo_interface__
                hull_poly = dict(type='Feature', properties=dict(id=1))
                hull_poly['geometry'] = geojs_geom
                #print(hull_poly)
                my_feature = Feature(geometry=geojs_geom)
                featCollection["features"].append(my_feature)
            print(featCollection)
            pl.plot(x,y,'o', color='#f16824')
        #
        #print(concave_hull.geoms)
    else:
        fig = pl.figure(figsize=(10,10))
        ax = fig.add_subplot(111)
        margin = .3
        x_min, y_min, x_max, y_max = concave_hull.bounds
        ax.set_xlim([x_min-margin, x_max+margin])
        ax.set_ylim([y_min-margin, y_max+margin])
        patch = PolygonPatch(concave_hull, fc='#999999', ec='#000000', fill=True, zorder=-1)
        ax.add_patch(patch)
        pl.plot(x, y, 'o', color='#f16824')
        geojs_geom = concave_hull.__geo_interface__
        my_feature = Feature(geometry=geojs_geom)
        featCollection["features"].append(my_feature)
    #
    file_object = open('e:\test\result\result.json', 'w')
    file_object.write(dumps(featCollection))
    file_object.close()
    subprocess.call(["ogr2ogr", "-t_srs", "EPSG:2385", "-f", "ESRI Shapefile", "e:\test\result\ConcaveEnvelope.shp", "e:\test\result\result.json"])
    pl.show()

    image

    image

  • 相关阅读:
    HDU 1203 【01背包/小数/概率DP】
    HDU 2955 【01背包/小数/概率DP】
    2018 计蒜之道 初赛 第三场
    判断一个多边形是顺时针还是逆时针的方法
    牛客网练习赛18 A 【数论/整数划分得到乘积最大/快速乘】
    Codeforces Round #311 (Div. 2)
    暑期训练狂刷系列——Hdu 3506 Largest Rectangle in a Histogram (单调栈)
    暑期训练狂刷系列——poj 3468 A Simple Problem with Integers (线段树+区间更新)
    暑期训练狂刷系列——Foj 1894 志愿者选拔 (单调队列)
    暑期训练狂刷系列——poj 3264 Balanced Lineup(线段树)
  • 原文地址:https://www.cnblogs.com/gispathfinder/p/5826236.html
Copyright © 2011-2022 走看看