Stinopys
Published in

Stinopys

Shortest distance to polygon on geoid surface

A fairly common geographic problem is finding the shortest distance from a reference point to some boundary line, or a point on a boundary line closest to a reference point. You can find many specific situations on your own that require a solution to this problem.

Searching the Internet for a solution to this problem led to an article of Maximilian HofmannCalculating distances from Points to Polygon Borders in Python — A Paris Example”. I liked the structure and presentation of the solution in the article, but I can not completely agree with author’s approach utilized finding the distance between the reference point and the middle of each segment of the polygon. Based on Maximilian’s approach, I used a more accurate method to solve the problem.

Content

Explanation of the method
The shortest distance between a point and a segment on the plane
Polygon on the plane
Spherical triangle
Spyder troubleshooting
Data preparation
Polygon
Coordinates of reference points
Description and call of the method
Displaying the result on the map
References
Ukrainian version

Explanation of the method

Now let’s return to our task.

The shortest distance between a point and a segment on the plane

The essence of finding the minimum distance between a point and a line segment on a plane can be explained without even using a picture. For this, let’s imagine a perpendicular CDdrawn from the reference point Cto the line, which is a continuation of the segment AB.

As a result, only 2 cases can occur:

  1. the point Dof intersection of the perpendicular and the line belongs to the segment AB — in this case, the minimum distance between the point C and the segment AB is equal to CD, i.e. the length of the perpendicular;
  2. the point Dof intersection of the perpendicular and the line does not belong to the segment AB — in this case, the minimum distance between the point Cand the segment ABis equal to CB(or CA, depending on which side of the segment the point is located D), i.e. the distance between the reference point C and the closest point from two extreme points of the segment — A and B;

Polygon on the plane

A polygon is in fact a connected sequence of line segments. Therefore, the task of finding the minimum distance between the reference point and the polygon is divided into separate, identical tasks of finding the minimum distance between the reference point and each line segment from the polygon. After finding such a minimum distance for each segment of the polygon, we choose the minimum value from the set of these distances.

Spherical triangle

In practice, the solution of geodesic problems is somewhat more complicated: the point and the polygon are not located on the plane, but on the surface of the geoid (Earth), which is sufficiently good described by the sphere. Then, instead of a flat triangle ABC, we have to consider a spherical triangle. Once before we have already dealt with such a triangle, in particular in the article “Mathematics of a spherical triangle”. In such a triangle, the sides are represented by arcs of a great circle of the sphere instead of straight segments. Curiously that main ratios of the spherical triangle are very similar to the ratios of the planar triangle, with the difference that the measure of the length of the first one is in angular units (degrees or radians), and not, say, in meters. Then, after the result for the spherical triangle is obtained in degrees, it can be converted into meters, knowing the radius of the sphere.

Spyder troubleshooting

Usually, the solution of any problem requires preliminary configuration of the programming environment. This case was not an exception due to the problem of improper operation of the Spyder 5.1.5 debugger with Python 3.9.6. To eliminate these problems, we had to install a virtual environment in Python with specific versions of the components, namely, by downgrading the version of Python from 3.9.6 to 3.8.10:

> conda create -n gis -c conda-forge python=3.8.10 spyder=5.1.5
folium geojson
> activate gis
> pip install PyGeodesy

Data preparation

Now we are ready for a practical solution to the problem.

As always, let’s start by importing the modules we need to solve the problem:

from pygeodesy import sphericalTrigonometry as st
import geojson
import folium
from folium.features import DivIcon

Test dataset for the construction of a polygon and set of reference points were created using Google Earth and then downloaded as kml files. For convenience, I converted data from kml to json format using one of the online converters, many of which can be found easily on the Internet. Downloading already converted test dataset in json format for the solution looks as follows:

## Read data
polygon_file = 'data/open-polygon.geojson'
points_file = 'data/points.geojson'

with open(polygon_file) as f:
data_poly = geojson.loads(f.read())
with open(points_file) as f:
data_pnts = geojson.loads(f.read())

Different methods, especially for geodesy purposes, use many formats to represent the same data. In our case, we will use two formats to represent our simple geodata:

  1. a list of lists of coordinates of individual points [Lat, Lon]and
  2. a list of the same LatLon coordinates needed to work with methods from the pygeodesy library .

To do this, we will perform another data conversion:

## Transforming data for 'nearestOn3' method
def transform_data(data):
# swap order of coords from [Long,Lat] to [Lat,Long]
pnts = [[p[1],p[0]] for p in
data.features[0].geometry.coordinates]
# prepare data type appropriate for 'nearestOn3' method
pnts_LL = [st.LatLon(p[0],p[1]) for p in pnts]
return pnts, pnts_LL

pnts, pnts_LL = transform_data(data_pnts)
poly, poly_LL = transform_data(data_poly)

⚠️ During these conversions, it is very important to control the correct order of coordinates (Lat, Lon)or (Lon, Lat)- in different applications and methods this order may be different!

Eventually I ended up using Google Earth to create arbitrary input test dataset with the following values:

Polygon

Coordinates of reference points

Description and call of the method

To find the nearest points and at the same time the shortest distance between the reference point and the polygon, we will use the nearestOn3 method from the PyGeodesy package . As mandatory input arguments, the method accepts the coordinates of the reference point pnt_LL and the list of polygon coordinates poly_LL. All coordinates are passed into the method under the type LatLon defined in the PyGeodesy package .

Given that we have defined several reference points through the list pnts_LL, let's immediately calculate in bulk the list of points closest nearest_pnts to these reference points:

nearest_pnts = [st.nearestOn3(pnt_LL, poly_LL) for pnt_LL in
pnts_LL]

Each element of the output list nearest_pnts is a tuple of three elements (closest, distance, angle), where:

  • closest– coordinates of the nearest point on the polygon poly_LL under the type LatLon or LatLon3Tuple(lat, lon, height);
  • distance– the distance between the nearest point closest on the polygon poly_LL and the reference point pnt_LL, which is converted to the same measurement units as the input argument radius, which is defined in meters by default radius = 6371008.77141;
  • angle is the beam angle from the reference point pnt_LL to the point closest in compass degrees.

For other purposes, particularly calculations and graphical display, we may need a list of nearest points nearest_pnts in a simpler format n_pnts, such as a list of lists of individual point coordinates [Lat, Lon], for that we perform a simple conversion:

n_pnts = [[np[0].lat,np[0].lon] for np in nearest_pnts]

Displaying the result on the map

Thanks to Maximilian Hofmann’s article, I saved a lot of effort in mastering the display of results. For that I reworked Maximilian’s code a bit.

The code snippet below for displaying the results is largely self-explanatory and consists of the following steps:

  1. determination of the central point poi of the map display window;
  2. display of the polygon data_poly on the map;
  3. display of reference points pnts and the corresponding nearest points n_pnts on the map;
  4. display of segments connecting the corresponding reference points pnts and the nearest points n_pnts on the map;
  5. saving the result in HTML format on disk.
#### DISPLAY RESULT

## Select central point for displaying the map
poi = data_pnts.features[0].geometry.coordinates[4]
poi = [poi[1], poi[0]] # in folium.Map the order is - Lat, Long

## Display polygon on the map
m = folium.Map(location=poi, zoom_start=9, tiles='openstreetmap')
folium.GeoJson(data_poly).add_to(m) # here the order is - Long, Lat

## Display reference points 'pnts' and nearest points 'n_pnts'
def show_pnts_on_map(points,map,color='red'):
fstr = '<div style="font-size: 16pt; color: ' + color + '">{}
</div>'
for i,_ in enumerate(points):
folium.CircleMarker(location=points[i], color=color,
radius=5, fill=color).add_to(map)
folium.map.Marker(location=points[i],
icon=DivIcon(icon_size=(150, 36),icon_anchor=(0, 0),
html=fstr.format(str(i+1)))).add_to(map)
return map

m = show_pnts_on_map(pnts,m,'red')
m = show_pnts_on_map(n_pnts,m,'blue')

## Display line segments connecting reference points and corresponding nearest points
for i,_ in enumerate(pnts):
folium.PolyLine(locations=[pnts[i],n_pnts[i]], color='red').add_to(m)
print(pnts[i][0], pnts[i][1], n_pnts[i][0], n_pnts[i][1],
nearest_pnts[i][1],nearest_pnts[i][2])

## Store the map on the disk
m.save('map_my_1.html')

You can view the result by simply clicking on the file map_my_1.htmlin the file explorer.

The result will look like this on the map:

In tabular form, the resulting dataset will look like this:

References

  1. Calculating distances from Points to Polygon Borders in Python — A Paris Example | by Maximilian Hofmann | Analytics Vidhya | Mar 19, 2020, Medium
  2. PyGeodesy — library of geodesic tools in Python
  3. Mathematics of the spherical triangle | Olegh Bondarenko | June 11, 2020, protw.io
  4. The code for finding the shortest distance from a point to a polygon on the surface of the geoid

Ukrainian version

--

--

Get the Medium app

A button that says 'Download on the App Store', and if clicked it will lead you to the iOS App store
A button that says 'Get it on, Google Play', and if clicked it will lead you to the Google Play store
Olegh Bondarenko

Olegh Bondarenko

Researcher, DSc, expert in Radiation Protection, Ecology, Air Quality Monitoring, Project Management, Data Science and other — orcid.org/0000-0001-8214-4654