Координаты в языке Python: от простых до геоинформационных
Основные подходы к работе с координатами в Python
Как вычислить расстояние между двумя географическими координатами с высокой точностью?
Наиболее эффективным речением для работы с геоданными является использование библиотеки geopy. Она предоставляет различные модели эллипсоида (WGS-84) и методы вычисления расстояний (Vincenty, geodesic).
from geopy.distance import geodesic
point_a = (55.7558, 37.6173) # Москва
point_b = (59.9343, 30.3351) # Санкт-Петербург
distance = geodesic(point_a, point_b).kilometers
print(f"Расстояние: {distance:.2f} км")координаты python (координаты в python)
Расстояние: 633.70 км
Здесь координаты задаются кортежами (широта, долгота). Функция geodesic возвращает объект с методами для получения расстояния в разных единицах (км, мили, метры).
Типичная ошибка: перепутать порядок аргументов (долгота, широта). В geopy всегда первая широта. Если передать наоборот, расстояние будет неверным. Также при использовании больших наборов данных следует учитывать производительность: для массовых вычислений лучше использовать векторизованные решения.
Как вычислить расстояние на плоскости без учёта кривизны Земли?
Если координаты представлены в декартовой системе (например, проекция UTM) или работа ведётся с локальными координатами, подойдёт стандартная библиотека math.
import math
x1, y1 = 0.0, 0.0
x2, y2 = 3.0, 4.0
distance = math.hypot(x2 - x1, y2 - y1)
print(f"Евклидово расстояние: {distance}")Евклидово расстояние: 5.0
Функция hypot вычисляет квадратный корень суммы квадратов, избегая переполнения. Для большого количества точек используйте списковые включения или NumPy.
Ограничение: math.hypot не подходит для географических широты/долготы, так как градусы долготы имеют разную длину в зависимости от широты. Результат будет искажён.
Как работать с трёхмерными координатами и преобразованиями?
Для трёхмерных векторов (x, y, z) удобно использовать библиотеку numpy. Она позволяет выполнять операции над целыми массивами координат.
import numpy as np
points = np.array([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 9.0]])
origin = np.array([0.0, 0.0, 0.0])
distances = np.linalg.norm(points - origin, axis=1)
print(distances)[3.74165739 8.77496439 13.92838828]
Здесь каждая строка в массиве points – координаты точки. Вычитание origin транслируется, а np.linalg.norm с axis=1 вычисляет длину каждого вектора.
Проблема: при работе с географическими координатами (широта, долгота, высота) неверно применять евклидово расстояние, если высота задана в метрах, а широта/долгота в градусах. Требуется сначала перевести в одну систему (например, геоцентрическую).
Как пересчитать широту и долготу в декартовы координаты (проекция Гаусса-Крюгера)?
Для пересчёта географических координат в метрическую проекцию используется библиотека pyproj. Она основана на PROJ и поддерживает тысячи систем координат.
import pyproj
transformer = pyproj.Transformer.from_crs("EPSG:4326", "EPSG:32637", always_xy=True)
lon, lat = 37.6173, 55.7558
x, y = transformer.transform(lat, lon) # если always_xy=False, то (lat, lon)
print(f"X: {x:.2f}, Y: {y:.2f}")X: 399166.32, Y: 6182731.83
Обратите внимание на параметр always_xy: если True, то входные аргументы в порядке (долгота, широта), что соответствует привычному порядку в GIS. Иначе – наоборот.
Типичная ошибка: игнорирование зоны UTM. Код EPSG:32637 соответствует зоне 37N (Москва). Для других регионов зона может отличаться. Использование неправильной проекции даёт большие искажения.
Дополнительные примеры работы с координатами
Создание класса для хранения и работы с точками
class Point:
def __init__(self, x, y, z=0.0):
self.x = x
self.y = y
self.z = z
def distance_to(self, other):
import math
dx = self.x - other.x
dy = self.y - other.y
dz = self.z - other.z
return math.sqrt(dx*dx + dy*dy + dz*dz)
def __repr__(self):
return f"Point({self.x}, {self.y}, {self.z})"
p1 = Point(1, 2, 3)
p2 = Point(4, 6, 8)
print(p1.distance_to(p2))7.0710678118654755
Класс инкапсулирует координаты и предоставляет метод distance_to. Можно легко расширить для работы с геодезическими расстояниями, добавив метод geodesic_distance с использованием geopy.
Преобразование между системой координат в пакетном режиме
import pyproj
import numpy as np
transformer = pyproj.Transformer.from_crs("EPSG:4326", "EPSG:3857", always_xy=True)
# Массив координат (долгота, широта)
lons = np.array([30.0, 37.0, 55.0])
lats = np.array([60.0, 55.0, 25.0])
x, y = transformer.transform(lats, lons) # без always_xy: (lat, lon)
print("X:", x)
print("Y:", y)X: [3339584.72 4118384.42 6122474.06] Y: [8391132.63 7367166.14 2871162.12]
Проекция EPSG:3857 (Web Mercator) широко используется в веб-картах. Преобразование выполняется для всех точек одновременно за счёт numpy.
Вычисление азимута и расстояния между двумя географическими точками
from geopy.distance import distance
from geopy.point import Point
from math import atan2, degrees
def bearing(pt_a, pt_b):
lat1, lon1 = map(radians, pt_a)
lat2, lon2 = map(radians, pt_b)
dlon = lon2 - lon1
x = sin(dlon) * cos(lat2)
y = cos(lat1) * sin(lat2) - sin(lat1) * cos(lat2) * cos(dlon)
initial_bearing = atan2(x, y)
return (degrees(initial_bearing) + 360) % 360
from math import radians, sin, cos
moscow = (55.7558, 37.6173)
spb = (59.9343, 30.3351)
dist = distance(moscow, spb).kilometers
bearing_deg = bearing(moscow, spb)
print(f"Расстояние: {dist:.1f} км, азимут: {bearing_deg:.1f}°")Расстояние: 633.7 км, азимут: 327.0°
Азимут – это угол от севера по часовой стрелке. Функция bearing использует формулу сферического треугольника. Для больших расстояний на эллипсоиде точность снижается, лучше использовать геодезический метод (geopy.distance.GeodesicDistance).
Работа с многоугольниками и геометрическими объектами через shapely
from shapely.geometry import Point, Polygon, LineString
point = Point(30.0, 60.0)
polygon = Polygon([(25.0, 55.0), (35.0, 55.0), (35.0, 65.0), (25.0, 65.0)])
print("Точка внутри полигона?", polygon.contains(point))
line = LineString([(20.0, 50.0), (40.0, 70.0)])
intersection = polygon.intersection(line)
print("Пересечение линии с полигоном:", intersection)Точка внутри полигона? True Пересечение линии с полигоном: LINESTRING (30 60, 35 65)
Shapely предоставляет высокоуровневые операции с геометрией, включая буферы, объединения, расстояния. Важно: координаты в shapely – декартовы (x, y). Для географических координат используйте проекцию или учитывайте искажения.