Algorytmy danych geoprzestrzennych

Jaki kształt ma Ziemia?

Geoida

Geoida

Geoida to model kształtu Ziemi, który przedstawia złożoną, nieregularną powierzchnię, na którą siła grawitacji działa wszędzie prostopadle do powierzchni.

Geoida nie ma matematycznie zdefiniowanego kształtu. Jej kształt wynika z różnic w ziemskim polu grawitacyjnym spowodowanymi zmianami w gęstości i rozkładzie masy Ziemi (np. występowanie gór i dolin).

Źródłem danych zazwyczaj są satelitarne pomiary grawitacyjne (np. misja GRACE).

Matematyczne reprezentacje

Elipsoida

Elipsoida to obiekt geometryczny, w którym wszystkie płaskie przekroje są elipsami lub okręgami. Wszystkie trzy osie mogą mieć różne długości.

\[ \frac{x^2}{a^2} + \frac{y^2}{a^2} + \frac{z^2}{b^2} = 1 \]

Elipsoida obrotowa (sferoida)

Elipsoida obrotowa to obiekt geometryczny utworzony przez obrót elipsy wokół jednej z jej osi. Jeśli elipsa zostanie obrócona wokół mniejszej osi, to w rezultacie powstanie elipsoida spłaszczona, która przypomina kształtem Ziemię.

Jest to szczególny przypadek elipsoidy, w której dwie z trzech osi są równe.

\[ \frac{x^2 + y^2}{a^2} + \frac{z^2}{c^2} = 1 \]

Sfera

Sfera to idealnie okrągły obiekt geometryczny, w którym każdy punkt na powierzchni jest w równej odległości od środka. Jest to podstawowe i ogólne przybliżenie kształtu Ziemi, ponieważ nie uwzględnia spłaszczenia na biegunach.

Jest to szczególny przypadek elipsoidy, w której wszystkie trzy osie mają równą długość.

\[ x^2 + y^2 + z^2 = r^2 \]

Obrazowanie Ziemi na mapie

Problem

Przedstawienie Ziemi na mapie obejmuje rzutowanie trójwymiarowego kształtu na dwuwymiarową powierzchnię. Ponieważ powierzchnia Ziemi jest zaokrąglona, proces ten wprowadza zniekształcenia związane ze kształtem, powierzchnią, odległością czy kierunkiem.

Odwzorowanie kartograficzne

Odwzorowanie kartograficzne to matematyczna transformacja szerokości i długości geograficznej z powierzchni sfery lub elipsoidy na płaszczyznę.

Podstawowe typy odwzorowań:

• Walcowe – Ziemia rzutowana jest na cylinder, który następnie jest rozwijany tworząc płaską mapę.
• Stożkowe – Ziemia rzutowana jest na stożek, który również jest rozwijany.
• Azymutalne (płaszczyznowe) – Ziemia rzutowana jest na płaską płaszczyznę, wyśrodkowaną w określonym punkcie.

Przykład

Przykład

Elipsa Tissota

Elipsa (wskaźnik) Tissota to matematyczna koncepcja używana w kartografii do opisu zniekształcenia kształtów i rozmiarów, które występuje podczas odwzorowania zakrzywionej powierzchni Ziemi na mapie.

Została opracowana przez francuskiego matematyka Nicolasa Auguste’a Tissota w XIX wieku.

Elipsa Tissota

W idealnym odwzorowaniu wskaźnik byłby okręgiem, wskazującym, że kształt i rozmiar obszaru są zachowane. Jednakże, w rzeczywistości wskaźniki te przyjmują kształt elips o różnych rozmiarach i orientacjach, w szczególności:

• jeśli oś wielka i mała elipsy nie są równe, to występuje zniekształcenie kształtu,
• jeśli rozmiar elipsy różni się od oryginalnego okręgu, to występuje zniekształcenie skali,
• jeśli osie elipsy nie są prostopadłe, to występuje zniekształcenie kątów (rotacja).

Elipsa Tissota

Wnioski

• Przedstawianie Ziemi na mapie to skomplikowany proces wymagający kompromisów.
• Każde odwzorowanie powoduje pewne zniekształcenia.
• Różne cele wymagają odmiennych podejść.
• Odwzorowanie jest zawsze przybliżeniem rzeczywistości.

Współrzędne geograficzne

Współrzędne geograficzne to zestaw dwóch liczb (szerokości i długości geograficznej) służący do określenia położenia punktu na powierzchni Ziemi. Mierzone są w stopniach, które są jednostkami miary kątowej.

Szerokość geograficzna (\(ϕ\)) określa, jak daleko na północ (N) lub południe (S) znajduje się lokalizacja od równika. Zakres od 0° na równiku do 90° na biegunach.

Długość geograficzna (\(λ\)) określa, jak daleko na wschód (E) lub zachód (W) znajduje się lokalizacja od południka zerowego. Zakres od 0° na południku zerowym do 180°.

Format zapisu

Zazwyczaj podawane są w jednym z dwóch formatów:

1. Stopnie (°), minuty (‘) i sekundy (’’). Jeden stopień podzielony jest na 60 minut, a jedna minuta na 60 sekund. Przykładowo: 52°17’34” N, 16°44’08” E.
2. Stopnie dziesiętne. Przykładowo: 52,2927° N, 16,7355° E.

Nie ma jednoznacznej reguły dotyczącej kolejności zapisu szerokości i długości geograficznej. W zależności od standardu mogą być różne!

Konwersja

Konwersja stopni, minut i sekund na stopnie dziesiętne:

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

Konwersja stopni dziesiętnych na stopnie, minuty i sekundy:

\[ \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 \]

Implementacja

decimal_to_dms = function(decimal_degrees, is_lat) {
  if (isTRUE(is_lat)) {
    if (decimal_degrees >= 0) direction = "N"
    else direction = "S" }
  else {
    if (decimal_degrees >= 0) direction = "E"
    else direction = "W" }
    
  decimal_degrees = abs(decimal_degrees)
  degrees = as.integer(decimal_degrees)
  minutes = as.integer((decimal_degrees - degrees) * 60)
  seconds = ((decimal_degrees - degrees) * 60 - minutes) * 60
  seconds = round(seconds)
  result = paste0(degrees, "°", minutes, "'", seconds, "''", " ", direction)
  return(result)
}

decimal_to_dms(52.2927, is_lat = TRUE)
#> "52°17'34'' N"
decimal_to_dms(16.7355, is_lat = FALSE)
#> "16°44'8'' E"

Implementacja

dms_to_decimal = function(dms_string) {
  dms_string = gsub("°", " ", dms_string)
  dms_string = gsub("'", " ", dms_string)
  parts = unlist(strsplit(dms_string, " "))
  parts = parts[parts != ""]
  
  degrees = as.numeric(parts[1])
  minutes = as.numeric(parts[2])
  seconds = as.numeric(parts[3])
  direction = parts[4]
  
  decimal = degrees + (minutes / 60) + (seconds / 3600)
  decimal = round(decimal, 4)
  if (direction %in% c("S", "W")) decimal = -decimal
  return(decimal)
}

dms_to_decimal("52°17'34'' N")
#> 52.2928
dms_to_decimal("16°44'08'' E")
#> 16.7356

Precyzja współrzędnych

Odległość

Rodzaje odległości

1. Odległość euklidesowa.
2. Odległość sferyczna.
3. Odległość elipsoidalna.
4. Odległość powierzchniowa.

Odległość euklidesowa

Odległość euklidesowa to miara odległości w linii prostej między dwoma punktami obliczana z następującego równania:

\[ 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.

Odległość euklidesowa

Odległość sferyczna

Sferyczne prawo cosinusów

\[ \Delta\sigma = \arccos\left(\sin(\phi_1)\sin(\phi_2) + \cos(\phi_1)\cos(\phi_2)\cos(\lambda_2 - \lambda_1)\right) \]

\[ d = r \cdot \Delta\sigma, \] gdzie:

\(\Delta\sigma\) to kąt środkowy między punktami na sferze,

\(r\) to średnia długość promienia Ziemi wynosząca 6371,009 km.

Wzór haversine

\[ a = \sin^2\left(\frac{\Delta\phi}{2}\right) + \cos (\phi_1) \cos (\phi_2) \sin^2\left(\frac{\Delta\lambda}{2}\right) \]

\[ c = 2 \cdot \operatorname{atan2}\left(\sqrt{a}, \sqrt{1-a}\right) \]

\[ d = r \cdot c, \]

gdzie:

\(\Delta\lambda = \lambda_2 - \lambda_1\) (różnica w odległości geograficznej),

\(\Delta\phi = \phi_2 - \phi_1\) (różnica w szerokości geograficznej),

\(r\) to średnia długość promienia Ziemi wynosząca 6371,009 km.

Odległość elipsoidalna

Odległość elipsoidalna to najkrótsza odległość między dwoma punktami na powierzchni elipsoidy. Jej obliczenie jest złożonym zadaniem i obejmuje rozwiązanie zestawu równań nieliniowych.

Najczęściej stosowaną metodą jest opracowana przez polsko-amerykańskiego geodetę Thaddeusa Vincenty’ego w 1975 r. Jednak, ograniczeniem mogą być problemy ze zbieżnością w przeciwległych punktach.

W przypadku zastosowań wymagających najwyższej precyzji stosuje się algorytm Karneya (2013 r.) uważany za najdokładniejszą metodę.

Odległość powierzchniowa

Odległość powierzchniowa

1. Wyodrębnij wartości wysokości terenu wzdłuż linii.
2. Oblicz skumulowaną odległość euklidesową w trzech wymiarach (X, Y, Z) między kolejnymi punktami.

points = extract(r, line, xy = TRUE)

total_length = 0
for (i in 2:nrow(points)) {
  p1 = c(points$X[i - 1], points$Y[i - 1], points$Z[i - 1])
  p2 = c(points$X[i], points$Y[i], points$Z[i])
  distance = sqrt(sum((p2 - p1)^2))
  total_length = total_length + distance
}

print(total_length)

Odległość powierzchniowa

Odległość planarna: 280 km

Odległość powierzchniowa: 283 km

Powierzchnia poligonów

Który z przedstawionych obszarów jest większy?

Powierzchnia poligonów

Powierzchnia poligonów

Każda kropka reprezentuje jeden kraj europejski.

Powierzchnia poligonów

Rzeczywisty rozmiar komórek

r = rast(xmin = 150000, ymin = 120000, xmax = 880000, ymax = 810000,
         resolution = 1000, crs = "EPSG:2180")
res(r)[1] * res(r)[2]
#> 1 000 000

Średnia powierzchnia komórek na sferze: 1 000 307 m\(^2\)

Porównanie modeli