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.

# Explanation of the method

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

1. 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;
2. 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 stimport geojsonimport foliumfrom 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 datapolygon_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' methoddef 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:

# 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 mappoi = data_pnts.features[0].geometry.coordinates[4]poi = [poi[1], poi[0]] # in folium.Map the order is - Lat, Long ​## Display polygon on the mapm = 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:

# Ukrainian version

--

--

## More from Stinopys

Search and research results are often lost. Although, no less often we have to return to them later. Graffiti (stinopys — in Ukrainian) is one of forms of abbreviated memorization and thus transformation of knowledge into valuable signposts.

## Olegh Bondarenko

42 Followers

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