Dokładniejszy sposób obliczania powierzchni rastrów

W mojej codziennej pracy ciągle jestem proszony o obliczanie obszarów globalnych zestawów danych rastrowych w projekcji geograficznej przy rozdzielczości 30 sekund kątowych. Te zestawy danych są zwykle wynikiem operacji Łączenia (typowym przykładem są klasy roślinności połączone z warstwą kraju). Aby to zrobić, nasza jednostka utworzyła zestaw danych rastrowych z obszarem każdego piksela w rzucie geograficznym z prędkością 30 sekund kątowych. Za pomocą tej siatki obszarów wykonywany jest status strefowy w celu zsumowania obszarów dla każdej klasy. Ponieważ nie jestem pewien, w jaki sposób utworzono tę siatkę obszaru, zawsze zastanawiałem się, czy to podejście jest bardziej dokładne po prostu przerzucenie rastra w rzucie o równej powierzchni (z prostych testów wyniki dwóch metod są podobne). Czy ktoś doświadczył podobnej sytuacji?

Istnieje stosunkowo prosta dokładna formuła dla obszaru dowolnego czworoboku kulistego ograniczonego równoległymi (liniami szerokości) i południkami (liniami długości). Można go wyprowadzić bezpośrednio przy użyciu podstawowych właściwości elipsy (osi głównej a i osi pomocniczej b ), która jest obracana wokół jej osi pomocniczej w celu wytworzenia elipsoidy. (Wyprowadzenie stanowi fajne ćwiczenie całkowe, ale uważam, że nie byłoby to interesujące na tej stronie).

Formułę upraszcza się, dzieląc obliczenia na podstawowe etapy.

Po pierwsze, odległość między wschodnią i zachodnią granicą - południki l0 i l1 - jest ułamkiem całego koła równym q = (l1 - 10) / 360 (gdy południki są mierzone w stopniach) lub 1 = ( l1 - 10) / (2 * pi) (gdy południki są mierzone w radianach). Znajdź obszar całego wycinka znajdujący się między podobieństwami f0 i f1 i pomnóż go przez q .

Po drugie, zastosujemy wzór na powierzchnię poziomego wycinka elipsoidy ograniczonej równikiem (dla f0 = 0) i równoległości na szerokości geograficznej f (= f1). Obszar plastra między dowolnymi dwoma szerokościami geograficznymi f0 i f1 (leżącymi na tej samej półkuli) będzie różnicą między większym i mniejszym obszarem.

Wreszcie, pod warunkiem, że model jest naprawdę elipsoidą (a nie kulą), obszar takiego przekroju między równikiem a równoleżnikiem na szerokości f jest określony przez

area(f) = pi * b^2 * (log(zp/zm) / (2*e) + sin(f) / (zp*zm))

gdzie ai bsą długościami odpowiednio osi głównej i pomocniczej elipsy generującej,

e = sqrt(1 - (b/a)^2)

jest jego ekscentryczność, i

zm = 1 - e*sin(f); zp = 1 + e*sin(f)

(Jest to o wiele prostsze niż obliczanie za pomocą geodezji, które i tak są jedynie przybliżeniem podobieństw. Zwróć uwagę na komentarz @cffk dotyczący sposobu obliczania log(zp/zm)w sposób pozwalający uniknąć utraty precyzji na małych szerokościach geograficznych).

area(f) to obszar nieprzezroczystego wycinka od równika do szerokości geograficznej f (na ilustracji około 30 stopni na północ. X i Y to geocentryczne kartezjańskie osie współrzędnych przedstawione w celach informacyjnych.

Dla elipsoidy WGS 84 użyj wartości stałych

a = 6 378 137 meters,  b =  6 356 752.3142 meters,

pociągający za sobą

e = 0.08181919084296

(Dla modelu sferycznej o = b , wzór się nieokreślone Trzeba się ograniczenie co e -.> 0 z góry, który następnie redukuje się do standardowego wzoru 2 * pi * a^2 * sin(f)).

Zgodnie z tymi wzorami czworokąt 30 na 30 w oparciu o równik ma powierzchnię 3077.2300079129 kilometrów kwadratowych, natomiast czworokąt 30 na 30 dotykający bieguna (który jest tak naprawdę tylko trójkątem) ma powierzchnię tylko 13,6086152 kwadratu kilometry.

Dla sprawdzenia, formuły zastosowane do wszystkich komórek siatki 720 na 360 pokrywającej powierzchnię ziemi dają całkowitą powierzchnię 4 * pi * (6371.0071809) ^ 2 kilometry kwadratowe, wskazując, że promień authaliczny ziemi powinien wynosić 6371,0071809 kilometrów. Różni się to od wartości Wikipedii tylko w ostatniej znaczącej liczbie (około jednej dziesiątej milimetra). (Myślę, że obliczenia Wikipedii są trochę wyłączone :-).

Jako dodatkowe kontrole wykorzystałem wersje tych wzorów do odtworzenia dodatków 4 i 5 w Lev M. Bugayevskiy i John P. Snyder, Map Projections: A Reference Manual (Taylor i Francis, 1995). W dodatku 4 pokazano długości łuku 30-metrowych odcinków południków i równoleżników, podanych do najbliższego metra. Kontrola wyników okazała się idealna. Następnie odtworzyłem tabelę z przyrostami 0,0005 ', zamiast przyrostów 0,5', i zintegrowałem numerycznie obszary czworoboku, jak oszacowano z tymi długościami łuku. Całkowity obszar elipsoidy został dokładnie odtworzony z dokładnością do ośmiu znaczących cyfr. Dodatek 5 pokazuje wartości area(f)dla f = 0, 1/2, 1, ..., 90 stopni pomnożone przez 1 / (2 * pi). Wartości te są podawane do najbliższego kilometra kwadratowego. Kontrola wzrokowa wartości bliskich 0, 45 i 90 stopni wykazała doskonałą zgodność.

Tę dokładną formułę można zastosować za pomocą algebry rastrowej, zaczynając od siatki podającej szerokości górne granice każdej komórki, a drugą podając szerokości geograficzne dolnych granic. Każda z nich jest zasadniczo siatką współrzędnych y. (W każdym przypadku może aby utworzyć sin(f)i zmi zpjako wyniki pośrednie.) Odjąć dwa wyniki, ma wartość bezwzględną, która, i mnożąc przez frakcji Q otrzymanego w pierwszym etapie (równe 0,5 / 360 = 1/720 na przykład dla komórki o szerokości 30 stóp). Będzie to siatka, której wartości zawierają dokładną wartośćobszary każdej komórki (do własnej dokładności liczbowej siatki). Upewnij się tylko, że wyrażasz szerokości geograficzne w postaci oczekiwanej przez funkcję sinusoidalną: wiele kalkulatorów rastrowych da ci współrzędne w stopniach, ale oczekujesz radianów dla ich funkcji trig!

Dla przypomnienia, tutaj są dokładne obszary komórek 30 na 30 na elipsoidzie WGS 84 od równika do bieguna, w odstępach 30 ', do 11 cyfr (ta sama liczba użyta dla mniejszego promienia b ):

3077.2300079,3077.0019391,3076.5458145,3075.8616605,3074.9495164,3073.8094348,3072.4414813,3070.8457347,3069.0222870,3066.9712434,3064.6927222,3062.1868550,3059.4537865,3056.4936748,3053.3066912,3049.8930202,3046.2528597,3042.3864209,3038.2939285,3033.9756204,3029.4317480,3024.6625762,3019.6683833,3014.4494612,3009.0061153,3003.3386648,2997.4474422,2991.3327939,2984.9950800,2978.4346744,2971.6519646,2964.6473522,2957.4212526,2949.9740951,2942.3063230,2934.4183938,2926.3107788,2917.9839636,2909.4384482,2900.6747464,2891.6933866,2882.4949115,2873.0798782,2863.4488581,2853.6024374,2843.5412166,2833.2658109,2822.7768503,2812.0749792,2801.1608571,2790.0351582,2778.6985716,2767.1518013,2755.3955665,2743.4306011,2731.2576543,2718.8774905,2706.2908892,2693.4986451,2680.5015685,2667.3004848,2653.8962347,2640.2896746,2626.4816763,2612.4731271,2598.2649300,2583.8580035,2569.2532818,2554.4517149,2539.4542684,2524.2619238,2508.8756783,2493.2965451,2477.5255533,2461.5637477,2445.4121891,2429.0719545,2412.5441367,2395.8298444,2378.9302026,2361.8463521,2344.5794500,2327.1306692,2309.5011988,2291.6922441,2273.7050264,2255.5407830,2237.2007674,2218.6862492,2199.9985139,2181.1388633,2162.1086151,2142.9091030,2123.5416769,2104.0077025,2084.3085615,2064.4456516,2044.4203864,2024.2341953,2003.8885234,1983.3848318,1962.7245972,1941.9093120,1920.9404843,1899.8196375,1878.5483108,1857.1280585,1835.5604507,1813.8470724,1791.9895239,1769.9894206,1747.8483931,1725.5680867,1703.1501618,1680.5962932,1657.9081707,1635.0874985,1612.1359952,1589.0553936,1565.8474409,1542.5138984,1519.0565410,1495.4771578,1471.7775513,1447.9595378,1424.0249466,1399.9756206,1375.8134157,1351.5402005,1327.1578567,1302.6682785,1278.0733724,1253.3750574,1228.5752643,1203.6759360,1178.6790272,1153.5865040,1128.4003439,1103.1225355,1077.7550785,1052.2999830,1026.7592702,1001.1349711,975.42912705,949.64378940,923.78101904,897.84288636,871.83147097,845.74886152,819.59715539,793.37845851,767.09488512,740.74855748,714.34160569,687.87616739,661.35438752,634.77841811,608.15041795,581.47255240,554.74699308,527.97591765,501.16150951,474.30595754,447.41145586,420.48020351,393.51440422,366.51626611,339.48800143,312.43182627,285.34996030,258.24462644,231.11805066,203.97246162,176.81009042,149.63317034,122.44393648,95.244625564,68.037475592,40.824725575,13.608615243

Wartości podano w kilometrach kwadratowych.

Jeśli chcesz przybliżyć te obszary lub po prostu lepiej zrozumieć ich zachowanie, formuła redukuje się do szeregu mocy według następującego wzoru:

area(f) = 2 * pi * b^2 * z * (1 + (4/3)y + (6/5)y^2 + (8/7)y^3 + ...)

gdzie

z = sin(f), y = (e*z)^2.

(Analogiczna formuła występuje w Bugayevskiy & Snyder, op. Cit. , Równanie (2.1).)

Ponieważ e ^ 2 jest tak małe (około 1/150 dla wszystkich elipsoidalnych modeli ziemi), a z leży między 0 a 1, y jest również małe. Tak więc terminy y ^ 2, y ^ 3, ... szybko się zmniejszają, dodając precyzję o kolejne dwa miejsca dziesiętne z każdym terminem. Gdybyśmy całkowicie zignorowali y , formuła byłaby wzorem pola kuli o promieniu b . Pozostałe terminy można rozumieć jako korygujące wybrzuszenie równikowe Ziemi.

Edytować

Pojawiły się pytania dotyczące porównania obliczeń odległości geodezyjnej z tymi dokładnymi wzorami. Metoda odległości geodezyjnej aproksymuje każdy czworokąt za pomocą geodezji, a nie równoległości, które łączą jego rogi poziomo, i stosuje wzór euklidesowy dla trapezu. W przypadku małych czworokątów, takich jak czworokąty 30 ', jest to nieco tendencyjne i ma względną dokładność między 6 a 10 części na milion. Oto wykres błędu dla WGS 84 (lub jakiejkolwiek uzasadnionej elipsoidy ziemskiej, jeśli o to chodzi):

Tak więc, jeśli (1) masz łatwy dostęp do obliczeń odległości geodezyjnej i (2) tolerujesz błąd na poziomie ppm, możesz rozważyć użycie tych obliczeń geodezyjnych i pomnożenie ich wyników przez 1,00000791 w celu skorygowania błędu systematycznego. W przypadku dwóch kolejnych miejsc dokładności dziesiętnej odejmij współczynnik korygujący pi / 2 * cos (2f) / 10 ^ 6: wynik będzie dokładny z dokładnością do 0,04 ppm.

Odpowiedź na pytanie radouxju zależy od kształtu piksela rzutowanego na elipsoidę. Jeśli układ współrzędnych rastra ma długość i szerokość geograficzną, piksel jest prostokątem linii loksodermy i można zastosować odpowiedź Whubera, lub, bardziej ogólnie, można użyć wzoru na wielokąt, którego krawędzie są liniami loksodromy. Jeśli układ współrzędnych jest rzutem konformalnym na dużą skalę (UTM, płaszczyzna stanu itp.), Dokładniejsze byłoby przybliżenie krawędzi za pomocą geodezji i zastosowanie wzoru na wielokąt geodezyjny. Wieloboki geodezyjne są prawdopodobnie najlepsze do ogólnego użytku, ponieważ w przeciwieństwie do wielokątów w loksodromie są „dobrze wychowane” blisko biegunów.

Implementacje formuł dla wielokątów geodezyjnych i loksodromów są dostępne w mojej bibliotece GeographicLib . Obszar geodezyjny jest dostępny w kilku językach; obszar linii rumbu to tylko C ++. Jest to wersja on-line (linia geodezyjna + Rhumb) dostępny tutaj . Dokładność tych obliczeń jest zwykle lepsza niż 0,1 metra kwadratowego.

Będziesz musiał osądzić na podstawie wiarygodnego / oficjalnego ... Wzory geodezyjne pochodzą z obszaru geodezyjnego (Danielsen, 1989, wymagana subskrypcja) oraz algorytmów geodezyjnych (Karney, 2013, otwarty dostęp). Wzory rumbu podano tutaj .

Natknąłem się na to pytanie, próbując określić wzór na obszar piksela WGS84. Podczas gdy odpowiedź @ Whuber zawiera te informacje, wciąż było trochę pracy, aby uzyskać formułę dla pola piksela w stopniach kwadratowych na danej szerokości geograficznej. Dołączyłem napisaną poniżej funkcję Pythona, która streszcza to w jednym wywołaniu. Chociaż nie odpowiada bezpośrednio na pytanie plakatu dotyczące obszaru CAŁEGO rastra (choć można by zsumować pola wszystkich pikseli), myślę, że nadal jest to użyteczna informacja dla kogoś, kto szuka podobnych obliczeń.

def area_of_pixel(pixel_size, center_lat):
    """Calculate m^2 area of a wgs84 square pixel.

    Adapted from: https://gis.stackexchange.com/a/127327/2397

    Parameters:
        pixel_size (float): length of side of pixel in degrees.
        center_lat (float): latitude of the center of the pixel. Note this
            value +/- half the `pixel-size` must not exceed 90/-90 degrees
            latitude or an invalid area will be calculated.

    Returns:
        Area of square pixel of side length `pixel_size` centered at
        `center_lat` in m^2.

    """
    a = 6378137  # meters
    b = 6356752.3142  # meters
    e = math.sqrt(1 - (b/a)**2)
    area_list = []
    for f in [center_lat+pixel_size/2, center_lat-pixel_size/2]:
        zm = 1 - e*math.sin(math.radians(f))
        zp = 1 + e*math.sin(math.radians(f))
        area_list.append(
            math.pi * b**2 * (
                math.log(zp/zm) / (2*e) +
                math.sin(math.radians(f)) / (zp*zm)))
    return pixel_size / 360. * (area_list[0] - area_list[1])