Want to share your content on python-bloggers? click here.
I’ve been spending time exploring Ningchuan Xiao’s excellent book, GIS Algorithms. This text offers a comprehensive dive into fundamental geospatial concepts including geometric algorithms, spatial indexing, and spatial analysis. The author thoroughly explains the theoretical foundations of each concept and provides practical Python implementations as examples.
In a series of upcoming posts, I’ll break down and share insights on these concepts one-by-one. In this post, we’ll focus on two straightforward geometric algorithms: calculating point-to-point and point-to-line distances. The accompanying code leverages the Shapely library to handle geometric primitives such as points and lines and to validate our custom implementations of each algorithm.
Point-to-Point Distance
We begin by computing the distance between two points in a rectilinear coordinate system. This can be accomplished using the Pythagorean theorem. For
In the next cell, point2point computes the distance between two points using the Pythagorean theorem:
from math import sqrt
from shapely.geometry import Point
def point2point(p0, p1):
"""
Return the distance between p1 and p2.
Parameters
----------
p0: shapely.geometry.point.Point
First point.
p1: shapely.geometry.point.Point
Second point.
Returns
-------
float
Straight-line distance between p0 and p1.
"""
x0: float = p0.x
y0: float = p0.y
x1: float = p1.x
y1: float = p1.y
return sqrt((x1 - x0)**2 + (y1 - y0)**2)
# Define two points.
p0 = Point(2., 3.)
p1 = Point(4., 7.)
# Compute distances using point2point and shapely.
p2p = point2point(p0, p1)
shp = p0.distance(p1)
print(f"Distance using point2point: {p2p:.5f}")
print(f"Distance using shapely : {shp:.5f}")
Distance using point2point: 4.47214 Distance using shapely : 4.47214
If it is necessary to compute distance between two points on a sphere separated by a large extent, using point2point will not be accurate. We can instead opt for the haversine formula, which calculates the great-circle distance between two points on the Earth’s surface. It is especially suited for applications where accuracy over large distances is critical. The haversine function in the next cell accepts two points, with x representing longitude and y latitude, and returns the great circle distance between them. To use the haversine formula, it is necessary to specify a value for the radius of the Earth. We use 6,371 km, a generally accepted value for the Earth’s mean radius (see here):
from math import asin, cos, sin, radians
def get_haversine(p0, p1):
"""
Calculate the great circle distance between two points on the earth
(specified in decimal degrees).
Parameters
----------
p0: shapely.geometry.point.Point
First point with x = longitude and y = latitude.
p1: shapely.geometry.point.Point
Second point with x = longitude and y = latitude.
Returns
-------
float
Great circle distance between p0 and p1 in km.
References
----------
http://www.movable-type.co.uk/scripts/latlong.html
"""
R: float = 6371.
lon0: float = p0.x
lat0: float = p0.y
lon1: float = p1.x
lat1: float = p1.y
# Convert lat and lon degrees to radians.
rlon0, rlat0, rlon1, rlat1 = map(radians, (lon0, lat0, lon1, lat1))
dlon = rlon1 - rlon0
dlat = rlat1 - rlat0
a = sin(dlat / 2)**2 + cos(rlat0) * cos(rlat1) * sin(dlon / 2)**2
c = 2 * asin(sqrt(a))
return R * c
The pyproj library can be used to validate our haversine distance calculation. p0 represents Chicago, Illinois and p1 Burlington, Vermont. pyproj.inv returns the distance in meters, so we divide by 1000 to covert to km:
import pyproj
# p0 = Chicago, IL, p1 = Burlington, VT.
p0 = Point(-87.6298, 41.8781)
p1 = Point(-73.2121, 44.4759)
# Compute distances using get_haversine and geopy.
get_haversine_dist = get_haversine(p0, p1)
geodesic = pyproj.Geod(ellps="WGS84")
_, _, pyproj_dist = geodesic.inv(*p0.coords[0], *p1.coords[0])
print(f"Distance (km) via get_haversine: {get_haversine_dist:,.0f}")
print(f"Distance (km) via pyproj : {pyproj_dist / 1000:,.0f}")
Distance (km) via get_haversine: 1,203 Distance (km) via pyproj : 1,205
The results are very close. We could get them to match exactly by tweaking the value used for the Earth’s radius.
Point-to-Line Distance
In Chapter 2 of GIS Algorithms, a general expression for the distance from a point to a line in rectilinear coordinates is presented. The distance
where
Let
resulting in
Squaring both sides and rearranging yields
After adding
and the R.H.S. as
Since the intersection point
In order to obtain the values for
where
Thus, a general form for the line can be expressed as
Multiplying through by
Referring back to our original expression for
Given a point
from shapely.geometry import LineString
def point2line(p, l):
"""
Compute distance between point p and line l.
Parameters
----------
p: shapely.geometry.point.Point
Point of interest.
l: shapely.geometry.linestring.LineString
Line defined by two end points.
Returns
-------
float
Distance from p to l.
"""
# Get endpoints from l.
p1, p2 = l.coords[:2]
x0: float = p.x
y0: float = p.y
x1: float = p1[0]
y1: float = p1[1]
x2: float = p2[0]
y2: float = p2[1]
dx = x1 - x2
dy = y1 - y2
a = dy
b = -dx
c = y1 * dx - x1 * dy
return abs(a * x0 + b * y0 + c) / sqrt(a**2 + b**2)
We create a line and point and use point2line to calculate the distance. The distance is also calculated using shapely to for verification:
p = Point(8., 3.)
l = LineString([[1., 1.], [5., 10.]])
point2line_dist = point2line(p, l)
shapely_dist = p.distance(l)
print(f"point2line distance: {point2line_dist:.5f}")
print(f"shapely distance : {shapely_dist:.5f}")
point2line distance: 5.58440 shapely distance : 5.58440
Let’s visualize p and l, along with the line perpendicular to l passing through p:
import matplotlib.pyplot as plt
# Determine slope and intercept of line perpendicular to l passing through p.
x0, y0 = p.coords[0]
(x1, y1), (x2, y2) = l.coords[:2]
# Get slope and intercept of l.
m = (y2 - y1) / (x2 - x1)
b = y1 - m * x1
# Get slope of line perpendicular to l.
mp = -1 / m
# Get intercept of line perpendicular to l that passes through (x0, y0).
bp = y0 - mp * x0
# Determine x and y where l and line perpendicular to l intersect.
# The two points that define the line perpendicular to l are (x0, y0)
# and (x_intersect, y_intersect).
x_intersect = (bp - b) / (m - mp)
y_intersect = m * x_intersect + b
# Group coordinates for plotting.
xx, yy = [x1, x2], [y1, y2]
xxp, yyp = [x0, x_intersect], [y0, y_intersect]
fig, ax = plt.subplots(1, 1, figsize=(7.5, 5.), tight_layout=True)
# Plot point p.
ax.plot(x0, y0, "o", markersize=8, color="red", label="original point")
# Plot line segment l.
ax.plot(xx, yy, linewidth=2.0, color="blue", label="original line")
# Plot line segment perpendicular to l.
ax.plot(xxp, yyp, linewidth=2.0, color="green", linestyle="--", label="perpendicular line")
# Plot intersection point.
ax.plot(x_intersect, y_intersect, "o", markersize=8, color="green", label="intersection point")
ax.annotate(
f"({x0:.2f}, {y0:.2f})", xy=(x0, y0), xycoords="data", ha="left", va="bottom",
textcoords="offset points", xytext=(5, 5), fontsize=10, rotation=0, weight="normal",
color="#000000"
)
ax.annotate(
"p", xy=(x0, y0), xycoords="data", ha="left", va="bottom", textcoords="offset points",
xytext=(-12, -12), fontsize=10, rotation=0, weight="normal", color="#000000"
)
ax.annotate(
f"({x_intersect:.2f}, {y_intersect:.2f})", xy=(x_intersect, y_intersect),
xycoords="data", ha="left", va="bottom", textcoords="offset points",
xytext=(5, 5), fontsize=10, rotation=0, weight="normal", color="#000000"
)
ax.axis('equal')
ax.tick_params(axis="x", which="major", direction='in', labelsize=7)
ax.tick_params(axis="y", which="major", direction='in', labelsize=7)
ax.xaxis.set_ticks_position("none")
ax.yaxis.set_ticks_position("none")
ax.legend(loc="upper right", fancybox=True, framealpha=1, fontsize="small")
plt.show()
Since we have the endpoints of the line perpendicular to l passing through p, we can show that point2point returns the same distance as point2line:
p1 = Point(x0, y0) p2 = Point(x_intersect, y_intersect) point2point(p1, p2)
5.584403908234904
Want to share your content on python-bloggers? click here.