# Algorithm for point to polygon distance calculations in shapely

I want to calculate the distance from a point to a polygon. If I run

``````point.distance(polygon)
``````

Will it:

1. Search the points that compose the exterior boundary of the polygon, select the closest point and return the distance to the original point
2. Search the vertices that compose the exterior boundary of the polygon, select that which has the closest node to the original point, then interpolate a new node at a position along the vertex that will represent the closest possible geometry in the polygon. Then calculate the distance from the original point to that node

I'm asking this because I'm trying to find the number of polygons that intersect a given radius of a point. I have been calculating Euclidean distance manually, searching the points in the exterior boundary of the polygon, but realise this might not be a representative distance (ie. accounting for the case where vertices that are within the radius while the nodes are not). I want my boundary method for calculating distances to match standard spatial libraries eg. shapely.distance

• It's open source so you can read the code to find out what it does exactly. Commented May 7 at 9:50
• True, but the distance calcs of shapely are calling a compiled instance of GEOS, where the process is a lot trickier to get in to - I figured there would be documentation about these operations? Commented May 7 at 9:54
• You can search through the GEOS code too, but in fact you want to be looking inside JTS to be sure. Commented May 7 at 9:56
• The GEOS documentation has some details but don't answer my question - they use the terminology 'distance between the closest points of two geometries'. To me this implies it uses the boundary points only (what I wrote as 1) but I am not sure if I'm interpreting correctly Commented May 7 at 9:58
• No, it finds the actual point Commented May 7 at 9:59

The second option you mentioned is the case, so not only points are taken in account, also any position on any line segment.

For illustration: both the distance of the green point and the red point to the polygon in the following image is 1, as shown in the code snippet below.

Regarding the precision, it it just regular math being used to calculate, no explicit grid. This results in quite precise measurements: in the little test/sample in the code snippet the distance from a point to the polygon is only rounded to 0 if the distance is smaller than 16 decimal places.

Code snippet:

``````import shapely
import shapely.plotting as plotter
import matplotlib.pyplot as plt

poly = shapely.box(1, 0, 5, 4)
point1 = shapely.Point(0, 0)
point2 = shapely.Point(0, 2)
point3 = shapely.Point(0.9999999999999999, 0)
point4 = shapely.Point(0.99999999999999999, 0)
print(f"{poly.distance(point1)=}")
print(f"{poly.distance(point2)=}")
print(f"{poly.distance(point3)=}")
print(f"{poly.distance(point4)=}")
# poly.distance(point1)=1.0
# poly.distance(point2)=1.0
# poly.distance(point3)=1.1102230246251565e-16
# poly.distance(point4)=0.0

fig, ax = plt.subplots()
plotter.plot_polygon(poly, ax=ax)
plotter.plot_points(point1, ax=ax, color="green")
plotter.plot_points(point2, ax=ax, color="red")
plt.show()
``````
• Thanks that's super helpful! I get that there is no explicit grid to calculate the distance, but if so how does it find the position? To put in context, my distance calc was interpolating points every x m along the line and then testing the distance against each one. This can't be right as the distance would depend on whether x was 0.01 m or 0.1 m, for example. Is there a formula/process that more efficiently finds the nearest point to calculate distance to along the line? (ideally just the equation, not a shapely function) Commented May 7 at 20:10
• I suppose it is just the perpendicular distance? More info on e.g. wikipedia. Commented May 7 at 20:20
• I wonder, why do you want to implement this yourself if there is an exisiting function in a performance-optimized library like GEOS? Commented May 7 at 20:36
• My idea is to subset the polygon dataset based on a KDTree (scipy spatial) distance from point query (using points interpolated along polygon boundaries) and then calculate precise distance using shapely/geos on only the features that had interpolated points within the radius. So far this is giving me much better performance than attribute subsets (eg. town the polygon is in) or STRtree nearest queries. I imagine the best performance would be using GEOS directly, but I'm not good enough with low-level languages to do this. Does that make sense? Very open to other ideas Commented May 8 at 22:05

The broader goal for this question was clarified in a comment on my previous answer. This answer offers a possible solution for this broader goal. Probably it would be better/clearer if this would be asked in a seperate question...

A typical way to search all geometries within distance of other geometries is to buffer the geometries with the maximum distance so, for points, they become ~circles of the size of the search distance. Then you can use a regular Rtree index (e.g. STRtree) to find all polygons within distance efficiently.

Typically the above will be sufficient, but be aware that the buffers are approximating curves by lines. So if you need really exact results for some reason, you can buffer slightly larger than the search distance, then filter again with the exact distance calculated.

Code sample:

``````import geopandas as gpd

from shapely import box, Point
import matplotlib.pyplot as plt

# Some test data
polys = gpd.GeoDataFrame(
data={
"poly_name": ["poly_left", "poly_right"],
"geometry": [box(0, 0, 2, 2), box(4, 0, 6, 2)],
}
)
points = gpd.GeoDataFrame(
{
"point_name": ["point1", "point2", "point3", "point4"],
"geometry": [Point(3, 0), Point(3, 1), Point(-0.5, 1), Point(-1, 3)],
},
)

max_distance = 1

# If you want to be really exact, you can make the buffer distance a bit large
# and do additional filtering based on the exact distance to compensate for
# buffer creating approximate circles rather than real circles.
buffer_distance = max_distance * 1
points_buffer = points.copy()
points_buffer["geometry"] = points_buffer.geometry.buffer(max_distance)

# Join the buffered points with the polygons. Typically this result is already
# fine. This uses an rtree under the hood, so will be fast.
within_distance = points_buffer.sjoin(polys, predicate="intersects")

# If it really needs to be exact and/or you want to know the distances,
# calculate the distances and/or do additional filtering here.
within_distance_points = gpd.GeoSeries(
within_distance[["point_name"]].join(points[["geometry"]]).geometry
)
within_distance_polys = gpd.GeoSeries(
within_distance[["index_right"]].set_index("index_right").join(polys).geometry
)
within_distance["distance"] = within_distance_points.distance(
within_distance_polys, align=False
)
within_distance = within_distance.loc[within_distance["distance"] <= max_distance]

# Print result
print(within_distance[["point_name", "poly_name", "distance"]])
#   point_name   poly_name  distance
# 0     point1  poly_right       1.0
# 1     point2   poly_left       1.0
# 1     point2  poly_right       1.0
# 2     point3   poly_left       0.5

# Plot input data
fig, ax = plt.subplots()
polys.plot(ax=ax, edgecolor="blue", facecolor="none")
points.plot(
ax=ax,
color="red",
facecolor="none",
)
points_buffer.plot(ax=ax, edgecolor="red", facecolor="none")
plt.show()
``````

Input data used for illustration: