# 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 Hofmann “Calculating 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 `CD`

drawn from the reference point `C`

to the line, which is a continuation of the segment `AB`

.

As a result, only 2 cases can occur:

- the point
`D`

of 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; - the point
`D`

of intersection of the perpendicular and the line**does not belong to the****segment**`AB`

— in this case, the minimum distance between the point`C`

and the segment`AB`

is 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:

- a list of lists of coordinates of individual points
`[Lat, Lon]`

and - 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:

- determination of the central point
`poi`

of the map display window; - display of the polygon
`data_poly`

on the map; - display of reference points
`pnts`

and the corresponding nearest points`n_pnts`

on the map; - display of segments connecting the corresponding reference points
`pnts`

and the nearest points`n_pnts`

on the map; - 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.html`

in the file explorer.

The result will look like this on the map:

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

# References

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