Stinopys
Published in

Stinopys

Найменша відстань до полігону на поверхні геоїду

Доволі поширеною географічною задачею є знаходження найменшої відстані від референтної точки до якоїсь межевої лінії, або найближчої точки з межевої лінії до референтної точки. Ви самостійно можете знайти безліч конкретних ситуацій, де вимагається вирішення цієї задачі.

Пошук в інтернеті рішення задачі вивів на статтю Maximilian HofmannCalculating distances from Points to Polygon Borders in Python”. Мені сподобалася структура і подання рішення в статті, але я не можу повністю погодитися з авторським підходом, що використовує визначення відстані між опорною точкою і серединою кожного відрізку полігону. Відштовхуючись від підходу Максиміліана, я використав точніший метод для вирішення поставленої задачі.

Зміст

Пояснення метода
Мінімальна відстань між точкою і відрізком на площині
Полігон на площині
Сферичний трикутник
Усунення проблеми зі Spyder
Підготовка даних
Полігон
Координати референтних точок
Опис і виклик метода
Відображення результату на мапі
Посилання
Англомовна версія статті

Пояснення метода

Тепер повернімося до нашої задачі.

Мінімальна відстань між точкою і відрізком на площині

Суть знаходження мінімальної відстані між точкою і відрізком прямої на площині можна пояснити, навіть не використовуючи малюнок. Для цього треба уявити перпендикуляр CD, проведений з референтної точки C до лінії, що є продовженням відрізку AB.

В результаті можуть виникнути лише 2 випадки:

  1. точка D перетину перпендикуляру і лінії належить відрізку AB – у цьому разі мінімальна відстань між точкою C і відрізком AB дорівнює CD, тобто довжині перпендикуляру;
  2. точка D перетину перпендикуляру і лінії не належить відрізку AB – у цьому разі мінімальна відстань між точкою C і відрізком AB дорівнює CB (або CA, в залежності від того, з якої сторони від відрізку розташована точка D), тобто відстань між референтною точкою C та найближчою точкою з двох крайніх точок відрізка — A і B.

Полігон на площині

Полігон представляє собою зв’язані між собою відрізки прямої. Тому задача знаходження мінімальної відстані між референтною точкою і полігоном розділяється на окремі однакові задачі пошуку мінімальної відстані між референтною точкою і кожним відрізком прямої зі складу полігону. Після знаходження для кожного відрізку полігону такої мінімальної відстані, ми вибираємо мінімальне значення з множини цих відстаней.

Сферичний трикутник

На практиці вирішення геодезничних задач дещо складніше: точка і полігон розташовані не на площині, а на поверхні геоїда (Землі), що доволі точно описується сферою. Тоді замість плаского трикутника ABC ми маємо розглядати сферичний трикутник. Ми вже колись розбирались з таким трикутником, зокрема у статті “Математика сферичного трикутника”. У такому трикутнику сторони представлені дугами великого кола сфери замість прямих відрізків. Цікаво, що основні співвідношення сферичного трикутника дуже схожі на співвідношення плаского трикутника з тією різницею, що міра довжини першого вимірюється у кутових одиницях (градусах чи радіанах), а не, скажімо, в метрах. Вже потім результат, отриманий для сферичного трикутника в градусах, можна перевести в метри, знаючи радіус сфери.

Усунення проблеми зі Spyder

Зазвичай вирішення будь-якої задачі вимагає попереднього налагодження середовища. Даний випадок не став виключенням через проблему неналежної роботи дебагера Spyder 5.1.5 з Python 3.9.6 . Для усунення цих неполадок прийшлося встановити віртуальний простір в Python з конкретними версіями компонент, а саме, понизивши версію Python з 3.9.6 до 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

Підготовка даних

Тепер ми готові до практичного вирішення задачі.

Як завжди почнімо із підключення потрібних нам для вирішення модулів:

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

Тестові дані для побудови полігону і декількох референтних точок були створені з допомогою Google Earth і скачані у вигляді kml-файлів. Для зручності я конвертував дані з kml у json формат з допомогою одного з онлайн-конверторів, безліч яких можна знайти в інтернеті. Завантаження вже конвертованих тестових даних у json форматі для вирішення виглядає наступним чином:

## 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())

Різні методи, особливо для цілей геодезії, використовують багато форматів для представлення тих самих даних. В нашому випадку ми використовуватимемо для представлення наших простих геоданих два формати:

  1. список списків координат окремих точок [Lat, Lon] і
  2. список тих самих координат типу LatLon, що необхідно для роботи з методами з бібліотеки pygeodesy.

Для цього проведемо ще одну конвертацію даних:

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

⚠️ Під час цих конвертацій дуже важливо контролювати правильний порядок координат (Lat, Lon) або (Lon, Lat) - у різних застосунках і методах цей порядок може бути різним!

В підсумку я створив з допомогою Google Earth довільні тестові вхідні дані, що мають такі значення:

Полігон

Координати референтних точок

Опис і виклик метода

Для знаходження найближчих точок і одночасно найкоротшої відстані між референтною точкою і полігоном використаємо метод nearestOn3 з пакету PyGeodesy . У якості обов’язкових вхідних аргументів метод приймає координати референтної точки pnt_LL і список координат полігону poly_LL. Всі координати передаються під типом LatLon, що визначений у пакеті PyGeodesy .

З огляду на те, що ми визначили декілька референтних точок через список pnts_LL, обрахуємо одразу список найближчих точок nearest_pnts до цих референтних точок:

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

Кожний елемент списку nearest_pnts представляє собою кортеж (tuple) з трьох елементів (closest, distance, angle), де:

  • closest – координати найближчої точки на полігоні poly_LL під типом LatLon або LatLon3Tuple(lat, lon, height);
  • distance – відстань між найближчою точкою closest на полігоні poly_LL і референтною точкою pnt_LL, що конвертована в ті самі одиниці вимірювання, що й вхідний аргумент radius, який за замовчанням визначений у метрах radius = 6371008.77141;
  • angle – кут променю від референтної точки pnt_LL до точки closest у градусах компасу.

Для інших цілей, зокрема, розрахунків і графічного відображення нам може знадобитись список найближчих точок nearest_pnts у простішому форматі n_pnts, як-от список списків координат окремих точок [Lat, Lon], для чого ми проведемо просте перетворення:

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

Відображення результату на мапі

Завдячуючи статті Maximilian Hofmann я зекономив багато зусиль на освоєнні відображення результатів і дещо переробив код Максиміліана.

Представлений нижче фрагмент коду відображення результатів є значною мірою самозрозумілим і складається з декількох наступних кроків:

  1. визначення центральної точки poi вікна відображення мапи;
  2. відображення полігону data_poly на мапі;
  3. відображення референтних точок pnts і відповідних найближчих точок n_pnts на мапі;
  4. відображення відрізків, що поєднують відповідні референтні точки pnts і найближчі точки n_pnts, на мапі;
  5. збереження результату у форматі HTML на диску.
#### 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')

Переглянути результат можна просто клацнувши файл map_my_1.html у файловому провідникові.

Результат виглядатиме на карті таким чином:

У табличному вигляді результуючи дані виглядатимуть таким чином:

Посилання

  1. Calculating distances from Points to Polygon Borders in Python — A Paris Example | by Maximilian Hofmann | Analytics Vidhya | Mar 19, 2020, Medium
  2. PyGeodesy — бібліотека геодезичних інструментів на Python
  3. Математика сферичного трикутника | Олег Бондаренко | 11 червня 2020, protw.io
  4. Код знаходження найменшої відстані від точки до полігону на поверхні геоїду

Англомовна версія статті

--

--

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.

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