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