Download as Notebook

#!pip install shapely
import shapely;
%matplotlib inline
Collecting shapely
 Downloading (1.5MB)
100% |████████████████████████████████| 1.5MB 11.0MB/s ta 0:00:01
Installing collected packages: shapely
Successfully installed shapely-1.6.4.post2
# A sample to learn with polygons, WKT, and relations
# (c) 2019 M. Werner

# Probably, you need to run
#   pip install shapely
# or
#   pip3 install shapely
# before using this script

from matplotlib import pyplot as plt;
from matplotlib import patches
from shapely import geometry
import numpy as np;

# If SNAP is true, then only integer coordinates in [0,10] are taken

# A helper function to plot a polygon. In python, you need to create and add a shape to an axis.
def plot_polygon(ax, poly, color):
    x, y = poly.exterior.coords.xy
    points = np.array([x, y], np.float32).T
    polygon_shape = patches.Polygon(points, linewidth=1, edgecolor='k', facecolor=color, alpha=0.5)

# Calculate the 9-IM matrix as a matrix.
def de9im(geom1, geom2):
    return(np.array([x for x in geom1.relate(geom2)]).reshape(3,3))

# Plotting the polygon (ignoring inner rings)
# Create a random polygon
dt = np.float32
if SNAP:
    dt = np.int8

A = geometry.Polygon (
   np.random.uniform(0,10,6).reshape(-1,2).astype(dt) # create 3x2
B = geometry.Polygon (
   np.random.uniform(0,10,6).reshape(-1,2).astype(dt) # create 3x2

# This now creates a plot window and extracts the axis object (needed for polygon drawing)
fig, ax = plt.subplots(1)
# Draw both polygons
plot_polygon(ax, A,'r')
plot_polygon(ax, B,'b')

# Print some relational information
print("Some Relations")
print("DE-9IM (String): %s" % A.relate(B))
print("WKT(A)        : %s" % (str(A)))
print("WKT(B)        : %s" % (str(B)))
print("A equals     B: %s" % A.equals  (B))
print("A contains   B: %s" % A.contains(B))
print("A crosses    B: %s" % A.crosses (B))
print("A disjoint   B: %s" % A.disjoint(B))
print("A intersects B: %s" % (A.intersects(B)))
print("A overlaps   B: %s" % (A.overlaps  (B)))
print("A touches    B: %s" % (A.touches   (B)))
print("A within     B: %s" % (A.within    (B)))

# zoom and show the plot

[['F' 'F' '2']
 ['F' 'F' '1']
 ['2' '1' '2']]
Some Relations
DE-9IM (String): FF2FF1212
WKT(A)        : POLYGON ((9 7, 4 2, 1 2, 9 7))
WKT(B)        : POLYGON ((7 9, 1 3, 0 6, 7 9))
A equals     B: False
A contains   B: False
A crosses    B: False
A disjoint   B: True
A intersects B: False
A overlaps   B: False
A touches    B: False
A within     B: False