93 lines
3.7 KiB
Python
93 lines
3.7 KiB
Python
import alphashape
|
|
from shapely.geometry import MultiPoint, mapping # Use mapping for GeoJSON conversion
|
|
import matplotlib.pyplot as plt
|
|
import json # To print GeoJSON nicely
|
|
|
|
# 1. Define a set of 2D points (e.g., forming a 'C' shape)
|
|
# Using simple integer coordinates for clarity
|
|
points = [
|
|
(0, 0), (0, 1), (0, 2), (0, 3), (0, 4), (0, 5), # Left vertical bar
|
|
(1, 5), (2, 5), (3, 5), # Top horizontal bar
|
|
(3, 4), # Small indent
|
|
(1, 0), (2, 0), (3, 0), # Bottom horizontal bar
|
|
# Add some interior points to make it more interesting
|
|
(1, 1), (1, 4), (2, 2.5)
|
|
]
|
|
|
|
print(f"Original points count: {len(points)}")
|
|
|
|
# 2. Calculate the alpha shape
|
|
# If alpha is not specified (or None), alphashape tries to find a 'good' value.
|
|
# You can experiment by setting alpha explicitly, e.g., alpha=2.0
|
|
try:
|
|
alpha_shape_geom = alphashape.alphashape(points) # Let alphashape determine alpha
|
|
# alpha_shape_geom = alphashape.alphashape(points, alpha=1) # Example with explicit alpha
|
|
|
|
print(f"\nAlpha shape calculated. Type: {alpha_shape_geom.geom_type}")
|
|
|
|
# Print as Well-Known Text (WKT)
|
|
print("Alpha Shape (WKT):")
|
|
print(alpha_shape_geom.wkt)
|
|
|
|
# Print as GeoJSON Geometry dictionary
|
|
print("\nAlpha Shape (GeoJSON Geometry):")
|
|
print(json.dumps(mapping(alpha_shape_geom), indent=2))
|
|
|
|
except Exception as e:
|
|
print(f"\nError calculating alpha shape: {e}")
|
|
alpha_shape_geom = None
|
|
|
|
# 3. Calculate the convex hull for comparison
|
|
# Use Shapely's MultiPoint for this
|
|
multi_point = MultiPoint(points)
|
|
convex_hull_geom = multi_point.convex_hull
|
|
|
|
print(f"\nConvex hull calculated. Type: {convex_hull_geom.geom_type}")
|
|
print("Convex Hull (WKT):")
|
|
print(convex_hull_geom.wkt)
|
|
|
|
|
|
# 4. Visualization using Matplotlib (Optional but helpful)
|
|
if alpha_shape_geom and convex_hull_geom:
|
|
fig, ax = plt.subplots(figsize=(8, 8))
|
|
|
|
# Plot the original points
|
|
x_coords, y_coords = zip(*points)
|
|
ax.scatter(x_coords, y_coords, color='blue', label='Original Points', s=50, zorder=3)
|
|
|
|
# Plot the convex hull
|
|
if convex_hull_geom.geom_type == 'Polygon':
|
|
hull_x, hull_y = convex_hull_geom.exterior.xy
|
|
ax.plot(hull_x, hull_y, color='red', linestyle='--', linewidth=2, label='Convex Hull', zorder=1)
|
|
elif convex_hull_geom.geom_type == 'LineString': # For collinear points
|
|
hull_x, hull_y = convex_hull_geom.xy
|
|
ax.plot(hull_x, hull_y, color='red', linestyle='--', linewidth=2, label='Convex Hull', zorder=1)
|
|
|
|
|
|
# Plot the alpha shape
|
|
# Handle both Polygon and MultiPolygon results from alphashape
|
|
def plot_polygon(ax, poly, color, label):
|
|
x, y = poly.exterior.xy
|
|
ax.plot(x, y, color=color, linewidth=2, solid_capstyle='round', label=label, zorder=2)
|
|
for interior in poly.interiors:
|
|
x_int, y_int = interior.xy
|
|
ax.plot(x_int, y_int, color=color, linewidth=1, linestyle=':', solid_capstyle='round')
|
|
|
|
if alpha_shape_geom.geom_type == 'Polygon':
|
|
plot_polygon(ax, alpha_shape_geom, color='green', label='Alpha Shape')
|
|
elif alpha_shape_geom.geom_type == 'MultiPolygon':
|
|
print("\nAlpha shape is a MultiPolygon, plotting all parts.")
|
|
for i, poly in enumerate(alpha_shape_geom.geoms):
|
|
plot_polygon(ax, poly, color='green', label=f'Alpha Shape Part {i+1}' if i==0 else None)
|
|
|
|
|
|
ax.set_title('Alpha Shape vs Convex Hull')
|
|
ax.set_xlabel('X coordinate')
|
|
ax.set_ylabel('Y coordinate')
|
|
ax.legend()
|
|
ax.set_aspect('equal', adjustable='box') # Ensure correct aspect ratio
|
|
plt.grid(True)
|
|
plt.show()
|
|
else:
|
|
print("\nSkipping visualization due to missing geometry.")
|