Algorytmy danych geoprzestrzennych
Obliczenia geodezyjne

Krzysztof Dyba

Współrzędne geograficzne

Współrzędne geograficzne to zestaw wartości liczbowych używany do określania położenia dowolnego punktu na powierzchni Ziemi. Dwoma podstawowymi składowymi współrzędnych geograficznych są:

  • szerokość geograficzna (latitude; \(\phi\)) – określa jak daleko na północ lub południe znajduje się punkt od równika Ziemi. Przyjmuje wartości z przedziału od 0° do 90° na północ (wartości dodatnie) i od 0° do 90° na południe (wartości ujemne) na biegunach.
  • długość geograficzna (longitude; \(\lambda\)) – określa jak daleko na wschód lub zachód znajduje się punkt od południka zerowego, który przebiega przez Greenwich w Anglii. Przyjmuje wartości z przedziału od 0° do 180° na wschód (wartości dodatnie) i od 0° do 180° na zachód (wartości ujemne).

Uwaga! Kolejność współrzędnych geograficznych może być różna w zależności od systemu, przyjętego standardu czy oprogramowania, np. norma ISO 6709 definiuje kolejność (szerokość, długość), podczas gdy standardy Open Geospatial Consortium określają ją odwrotnie (patrz GeoJSON czy Well-Known Text).

Współrzędne geograficzne są zazwyczaj wyrażane w jednym z następujących formatów:

  • stopnie, minuty i sekundy (degrees, minutes, seconds; DMS) – stopnie reprezentują miarę kątową. Każdy stopień podzielony jest na 60 minut, a każda minuta na 60 sekund. Przykładowo: 52°17’34” N, 16°44’08” E.
  • stopnie dziesiętne (decimal degrees; DD) – współrzędne wyrażone są jako wartości dziesiętne stopni. Przykładowo: 52,2927° N, 16,7355° E.
  • współrzędne planarne – współrzędne wyrażone jako pomiary liniowe w układzie planarnym. Przykładowo: 345613,2393 m, 494269,4463 m (układ PUWG 1992).

Konwersja

Wartość DMS konwertowana jest na stopnie dziesiętne przy użyciu następującego wzoru:

\[ \text{Stopnie dziesiętne} = \text{Stopnie} + \frac{\text{Minuty}}{60} + \frac{\text{Sekundy}}{3600} \]

def dms_to_decimal(dms_str):
    
    # zamień znaki na spacje
    dms_str = dms_str.replace("°", " ").replace("'", " ").replace("\"", " ")
    # podziel tekst na części według spacji
    parts = dms_str.split()
    
    degrees = int(parts[0])
    minutes = int(parts[1])
    seconds = int(parts[2])
    direction = parts[3]
    
    decimal = degrees + (minutes / 60) + (seconds / 3600)
    decimal = round(decimal, 4)
    
    # sprawdź kierunek i nadaj znak
    if direction in ["S", "W"]:
        decimal = -decimal
    
    # zwracany jest typ liczbowy
    return decimal
print(dms_to_decimal("52°17'34\" N"), "°", sep = "")
print(dms_to_decimal("16°44'08\" E"), "°", sep = "")
52.2928°
16.7356°

Natomiast operacja odwrotna, czyli konwersja stopni dziesiętnych na wartość DMS, wymaga zastosowanie poniższych wzorów:

\[ \text{Stopnie dziesiętne} = |\text{Stopnie dziesiętne}| \]

\[ \text{Stopnie} = \text{int}(\text{Stopnie dziesiętne}) \]

\[ \text{Minuty} = \text{int}\left((\text{Stopnie dziesiętne} - \text{Stopnie}) \times 60\right) \]

\[ \text{Sekundy} = \left((\text{Stopnie dziesiętne} - \text{Stopnie}) \times 60 - \text{Minuty}\right) \times 60 \]

Obliczenia

Odległość

W tej sekcji porównamy różne metody obliczania odległości między punktami uwzględniając:

  • odległość euklidesową (układ kartezjański),
  • odległość sferyczną,
  • odległość elipsoidalną.

W tym celu wykorzystamy dwa punkty reprezentujące lokalizację Warszawy (21,0122°, 52,2296°) oraz Rzymu (12,5113°, 41,8919°).

Euklidesowa

Odległość euklidesowa w przestrzeni dwuwymiarowej definiowana jest jako:

\[ d=\sqrt{(\Delta x)^2+(\Delta y)^2}=\sqrt{(x_2-x_1)^2+(y_2-y_1)^2} \]

Wynika ona bezpośrednio z twierdzenia Pitagorasa o trójkątach prostokątnych:

\[ a^2 + b^2 = c^2, \]

gdzie \(a\) i \(b\) to długość przyprostokątnych, a \(c\) to długość przeciwprostokątnej.

Procedura wymaga zastosowania transformacji współrzędnych z układu elipsoidalnego do kartezjańskiego, aby uzyskać sensowne wyniki. Należy wykorzystać klasę QgsCoordinateTransform i na przykład odwzorowanie sferyczne Merkatora (EPSG:3857). Po przeprowadzonej transformacji, można obliczyć odległość między punktami używając metodę distance(). Zauważ, że to odwzorowanie znacząco zniekształca rozmiar i odległość, zatem bardziej prawidłowym podejściem byłoby zastosowanie odwzorowania równoodległościowego.

from qgis.core import QgsPointXY, QgsCoordinateReferenceSystem
from qgis.core import QgsCoordinateTransform, QgsCoordinateTransformContext

pt1 = (21.0122, 52.2296)
pt2 = (12.5113, 41.8919)

pt1 = QgsPointXY(pt1[0], pt1[1])
pt2 = QgsPointXY(pt2[0], pt2[1])

crs_4326 = QgsCoordinateReferenceSystem("EPSG:4326")
crs_3857 = QgsCoordinateReferenceSystem("EPSG:3857")
transform_context = QgsCoordinateTransformContext()
transform = QgsCoordinateTransform(crs_4326, crs_3857, transform_context)

pt1_3857 = transform.transform(pt1)
pt2_3857 = transform.transform(pt2)

distance = pt1_3857.distance(pt2_3857)
distance = round(distance / 1000, 2) # wynik w km

print("Odległość między Warszawą a Rzymem:", distance, "km")
Odległość między Warszawą a Rzymem: 1942.97 km

Możemy również samodzielnie zaimplementować metodę distance() dla obiektów typu QgsPointXY w następujący sposób:

def euclidean_distance(pt1, pt2):
    import math
        
    # wyodrębnienie współrzędnych z obiektu QgsPoint
    x1, y1 = pt1.x(), pt1.y()
    x2, y2 = pt2.x(), pt2.y()
        
    distance = math.sqrt( (x2 - x1)**2 + (y2 - y1)**2 )
    distance = round(distance / 1000, 2) # wynik w km
    return distance

euclidean_distance(pt1_3857, pt2_3857)
1942.97

Sferyczna

Drugie podejście opiera się na wykorzystaniu trygonometrii na sferze, która przybliża powierzchnię Ziemi. W tym przypadku pomijamy transformację układu współrzędnych. Jednakże, do obliczeń zamiast stopni należy wykorzystać radiany. Ten sposób oferuje błąd na poziomie 0,5%.

Wzór jest następujący:

\[ \Delta{x} = \cos(\phi_2)\cos(\lambda_2) - \cos(\phi_1)\cos(\lambda_1) \]

\[ \Delta{y} = \cos(\phi_2)\sin(\lambda_2) - \cos(\phi_1)\sin(\lambda_1) \]

\[ \Delta{z} = \sin(\phi_2) - \sin(\phi_1) \]

\[ d_\textrm{t} = r \sqrt{(\Delta{x})^2 + (\Delta{y})^2 + (\Delta{z})^2} \]

\[ d = 2 r \arcsin \frac{d_\textrm{t}}{2 r}, \]

gdzie \(r\) to średni promień Ziemi równy 6 371 009 m.

Średni promień Ziemi można wyliczyć z poniższego równania:

\[ r = \frac{2a+b}{3}, \]

gdzie \(a\) to promień równikowy wynoszący 6 378 137 m, a \(b\) to promień biegunowy wynoszący 6 356 752 m.

# https://en.wikipedia.org/wiki/Geographical_distance#Spherical-surface_formulae

def spherical_distance(pt1, pt2):
    from math import sin, cos, asin, sqrt, pi
    
    r = 6371009 # średni promień Ziemi w metrach
        
    coords = [pt1.x(), pt1.y(), pt2.x(), pt2.y()]
    # konwersja stopni na radiany
    coords = [coord * pi / 180 for coord in coords]
    x1, y1, x2, y2 = coords
    
    delta_x = cos(y2) * cos(x2) - cos(y1) * cos(x1)
    delta_y = cos(y2) * sin(x2) - cos(y1) * sin(x1)
    delta_z = sin(y2) - sin(y1)
    d_t = r * sqrt( delta_x**2 + delta_y**2 + delta_z**2 )
    distance = 2 * r * asin(d_t / (2 * r))
    distance = round(distance / 1000, 2)
    
    return distance
distance = spherical_distance(pt1, pt2)
print("Odległość między Warszawą a Rzymem:", distance, "km")
Odległość między Warszawą a Rzymem: 1315.51 km

Elipsoidalna

Precyzyjniejszy wynik można uzyskać stosując obliczenia na elipsoidzie zamiast sfery (wzór Vincenty’ego). Służy do tego klasa QgsDistanceArea oraz metoda measureLine(). Należy również wybrać elipsoidę używając metody setEllipsoid(). Standardowo do pomiarów globalnych wykorzystuje się elipsoidę WGS84.

from qgis.core import QgsDistanceArea

distance_area = QgsDistanceArea()
distance_area.setEllipsoid("WGS84")

distance = distance_area.measureLine(pt1, pt2)
distance = round(distance / 1000, 2) # wynik w km

print("Odległość między Warszawą a Rzymem:", distance, "km")
Odległość między Warszawą a Rzymem: 1316.2 km

Powierzchnia

Sprawdźmy teraz jakie różnice zachodzą w przypadku pomiarów powierzchni w układzie planarnym (EPSG:2180) oraz na elipsoidzie. Jako dane wejściowe wykorzystajmy przybliżoną obwiednię Polski.

Kartezjańska

Procedura wygląda identycznie jak w przypadku obliczania odległości między punktami z tą różnicą, iż należy wykorzystać metodę area().

from qgis.core import QgsGeometry

bbox = "POLYGON((14 49, 14 54, 24 54, 24 49, 14 49))"
geom = QgsGeometry.fromWkt(bbox)

crs_4326 = QgsCoordinateReferenceSystem("EPSG:4326")
crs_2180 = QgsCoordinateReferenceSystem("EPSG:2180")
transform_context = QgsCoordinateTransformContext()
transform = QgsCoordinateTransform(crs_4326, crs_2180, transform_context)

Uwaga! Transformacja geometrii jest dokonana bezpośrednio na obiekcie klasy QgsGeometry.

geom.transform(transform) # EPSG:2180
print(geom.asWkt(precision = 1))
Polygon ((134461.7 137878.5, 172479.3 693299.8, 827520.7 693299.8, 865538.3 137878.5, 134461.7 137878.5))
area = geom.area()
area = round(area / 1000**2, 2) # wynik w km^2

print("Powierzchnia poligonu wynosi:", area, "km^2")
Powierzchnia poligonu wynosi: 384939.74 km^2

Elipsoidalna

Do obliczenia powierzchni elipsoidalnej należy użyć metody measureArea(), tak jak zademonstrowano to w poprzedniej sekcji na przykładzie odległości elipsoidalnej.

geom = QgsGeometry.fromWkt(bbox) # EPSG:4326

distance_area = QgsDistanceArea()
distance_area.setEllipsoid("WGS84")

area = distance_area.measureArea(geom)
area = round(area / 1000**2, 2)

print("Powierzchnia poligonu wynosi:", area, "km^2")
Powierzchnia poligonu wynosi: 385343.21 km^2

Zadania

  1. Zaimplementuj funkcję decimal_to_dms() umożliwiającą konwersję stopni dziesiętnych na wartość DMS. W funkcji uwzględnij argument logiczny is_lat, który pozwoli określić czy współrzędna wejściowa reprezentuje szerokość czy długość geograficzną i zwróci wynik z odpowiednim kierunkiem geograficznym.
  2. Zaimplementuj wzór haversine [9] jako funkcję i porównaj z innymi otrzymanymi wynikami na przykładzie odległości pomiędzy Warszawą a Rzymem. Współrzędne geograficzne wyrażone w stopniach muszą zostać zamienione na radiany.
  3. Zaimplementuj funkcję surface_distance(), która obliczy odległość powierzchniową między punktami w przestrzeni trójwymiarowej. Następnie wykorzystaj ją do obliczenia całkowitej długości profilu terenu zapisanego w pliku profil.csv w układzie EPSG:2180. Do wczytania pliku csv można wykorzystać moduł csv lub klasę QgsVectorLayer w QGIS.