生成凹包

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()