Używasz SRTM Global DEM do obliczania nachylenia?

17

Pobrałem SRTM GDEM (rozdzielczość ~ 90 km).

Korzystam z ArcGIS 10.

Próbowałem użyć analityka przestrzennego do obliczenia nachylenia.

Nie mogę jednak obliczyć nachylenia.

Wartości wyjściowe mają tylko dwa zakresy 0 i 0,1–90.

Nie jestem pewien, na czym polega problem?

arcgis-desktop arcgis-10.0 dem srtm użytkownik2543
źródło
To zależy od tego, gdzie analizujesz na świecie. Istnieją różne prognozy dla każdej lokalizacji. Gdzie egzaminujesz
djq
5
Rozdzielczość wynosi w rzeczywistości ~ 90 m, a nie ~ 90 km.
Akheloes,
Tylko komentarz, jeśli jesteś na utrzymaniu pulpitu, możesz zalogować się do ArcGIS Online i korzystać z ich usług elewacji (bez potrzeby rozszerzenia NA). Warstwa nachylenia może być używana jako warstwa odniesienia. W Australii mamy 1 sekundowe dane SRTM (~ 30m res) blogs.esri.com/esri/arcgis/2014/07/11/…
Simon

Odpowiedzi:

29

Wydaje się, że to dobre miejsce na opisanie prostego, szybkiego i więcej niż dość dokładnego sposobu obliczania stoków dla globalnie rozbudowanego DEM .

Zasady

Przypomnijmy, że nachylenie powierzchni w punkcie jest zasadniczo największym stosunkiem „wznoszenia” do „biegu” spotykanym we wszystkich możliwych łożyskach z tego punktu. Problem polega na tym, że gdy rzut ma zniekształcenie skali, wartości „przebiegu” zostaną niepoprawnie obliczone. Co gorsza, gdy zniekształcenie skali zmienia się wraz z łożyskiem - co ma miejsce w przypadku wszystkich niezgodnych rzutów - to, jak nachylenie zmienia się w zależności od łożyska, zostanie niepoprawnie oszacowane, uniemożliwiając dokładną identyfikację maksymalnego współczynnika wzrostu: przebiegu (i przekrzywienia obliczenie aspektu).

Możemy to rozwiązać za pomocą rzutu konformalnego, aby upewnić się, że zniekształcenie skali nie zmienia się w zależności od namiaru, a następnie korygując oszacowania nachylenia, aby uwzględnić zniekształcenie skali (które zmienia się od punktu do punktu na mapie). Sztuką jest użycie globalnej projekcji konformalnej, która pozwala na proste wyrażenie zniekształcenia skali.

Rzut Mercatora pasuje do rachunku: zakładając, że skala jest równa na równiku, jego zniekształcenie jest równe siecznej szerokości geograficznej. Oznacza to, że odległości na mapie wydają się pomnożone przez siecznego. Powoduje to, że wszelkie obliczenia nachylenia obliczają wzrost: (sec (f) * run) (co jest stosunkiem), gdzie f jest szerokością geograficzną. Aby to naprawić, musimy pomnożyć obliczone nachylenia przez sekundę (f); lub równoważnie podziel je przez cos (f). To daje nam prosty przepis:

Oblicz nachylenie (jako wzrost: bieg lub procent) za pomocą rzutu Mercatora, a następnie podziel wynik przez cosinus szerokości geograficznej.

Przepływ pracy

Aby to zrobić z siatką podaną w stopniach dziesiętnych (np. SRTM DEM), wykonaj następujące czynności:

  1. Utwórz siatkę szerokości geograficznej. (To tylko siatka współrzędnych y).

  2. Oblicz jego cosinus.

  3. Rzutuj zarówno DEM, jak i cosinus szerokości geograficznej, używając rzutu Mercatora, w którym skala jest prawdziwa na równiku.

  4. W razie potrzeby przekonwertuj jednostki elewacji, aby zgadzały się z jednostkami rzutowanych współrzędnych (zwykle metrów).

  5. Oblicz nachylenie rzutowanego DEM jako czyste nachylenie lub procent ( nie jako kąt).

  6. Podziel to nachylenie przez rzutowaną siatkę cosinus (szerokość geograficzna).

  7. W razie potrzeby ponownie skieruj siatkę nachylenia na dowolny inny układ współrzędnych w celu dalszej analizy lub mapowania.

Błędy w obliczeniach nachylenia będą wynosić do 0,3% (ponieważ ta procedura wykorzystuje kulisty model ziemi zamiast modelu elipsoidalnego, który jest spłaszczony o 0,3%). Ten błąd jest znacznie mniejszy niż inne błędy, które przechodzą do obliczeń nachylenia i dlatego można go pominąć.

W pełni globalne obliczenia

Projekcja Mercator nie obsługuje żadnego bieguna. Do pracy w regionach polarnych należy rozważyć użycie polarnej projekcji stereograficznej z rzeczywistą skalą na biegunie. Zniekształcenie skali wynosi 2 / (1 + sin (f)). Użyj tego wyrażenia zamiast sec (f) w przepływie pracy. W szczególności, zamiast obliczać siatkę cosinus (szerokość geograficzna), oblicz siatkę, której wartości to (1 + sin (szerokość)) / 2 ( edytuj : użyj -latitude dla Bieguna Południowego, jak omówiono w komentarzach). Następnie postępuj dokładnie tak jak poprzednio.

Aby uzyskać kompletne globalne rozwiązanie, należy rozważyć podzielenie siatki naziemnej na trzy części - jedną wokół każdego bieguna i drugą wokół równika - wykonując obliczenia nachylenia osobno dla każdej części za pomocą odpowiedniego rzutu i mozaikując wyniki. Rozsądnym miejscem do podzielenia kuli ziemskiej są kręgi szerokości geograficznej na szerokości 2 * ArcTan (1/3), czyli około 37 stopni, ponieważ na tych szerokościach geograficznych współczynniki korekcji Mercator i stereograficzne są sobie równe (mają wspólną wartość 5/4) i byłoby miło zminimalizować rozmiary dokonanych poprawek. Aby sprawdzić obliczenia, siatki powinny być w bardzo ścisłej zgodności, w których się pokrywają (niewielkie ilości niedokładności zmiennoprzecinkowej i różnice wynikające z ponownego próbkowania rzutowanych siatek powinny być jedynymi źródłami rozbieżności).

Bibliografia

John P. Snyder, Projekcje map - Podręcznik roboczy . USGS Professional Paper 1395, 1987.

Whuber
źródło
7
Znajduję się w pozycji, jak to często robię, by jeszcze raz podziękować whuberowi za opisanie rozwiązania, a także za przedstawienie uzasadnienia, które go buduje. Mój kapelusz jest od pana.
matt wilkie
Dzięki @matt. Nie chciałem wcześniej sugerować, że twoja (teraz usunięta) odpowiedź powinna zostać odwołana: w rzeczywistości głosowałem za nią, ponieważ udostępniłeś link do interesującego odniesienia USGS, które może być przydatne dla wielu czytelników. (Mój komentarz był krytyczny tylko dla drugiego fragmentu tego artykułu, a nie samego artykułu.)
whuber
ach Dziękuję za wyjaśnienie. Odtworzyłem odpowiedź, ufając, że ludzie mają teraz przed sobą wystarczającą ilość informacji, aby dokonać świadomego wyboru :)
Matt Wilkie
2
Pochodząc z francuskiego tła, zajęło mi trochę czasu, aby przetłumaczyć niezbędną terminologię, aby lepiej zrozumieć tę świetną odpowiedź, więc pomyślałem, że upuszczenie tego linku jest dobrą pomocą dla początkujących, takich jak ja: webhelp.esri.com/arcgisdesktop/9.2/...
Akheloes
To świetne podejście i już wykorzystałem twoje rozwiązanie do wygenerowania globalnego rastra nachyleń. Jedna wskazówka z praktycznego doświadczenia: Ponieważ wartości szerokości geograficznej na południe od równika są ujemne, musisz użyć bezwzględnej wartości szerokości geograficznej w następującym równaniu: (1 + sin (szerokość geograficzna)) / 2
Saleika
18

Oryginalna odpowiedź

Chyba poziome jednostki dla twojego rastra są w stopniach lub sekundach kątowych. Musisz przerzucić ten raster na rzut przestrzenny, w którym twoje jednostki poziome i pionowe są takie same (tj. Jeśli jednostki pionowe są w metrach, sugeruję użycie UTM, który ma poziome jednostki metrów).

Aby zmienić projekt rastra za pomocą ArcCatalog / ArcGIS, zajrzyj do:

ArcToolbox> Narzędzia do zarządzania danymi> Projekcje i transformacje> Raster> Project Raster

Wybierz rzutowane odniesienie przestrzenne obejmujące twój region zainteresowania, np. Wypróbuj strefę UTM. Istnieje wiele innych opcji, które najlepiej udokumentować w instrukcji . Uwaga: nie można utworzyć zbioru danych nachylenia dla całej Ziemi (jeśli to właśnie próbujesz zrobić).

Lepsza odpowiedź, używając GDAL ze skalą

Teraz, gdy dane SRTM są globalnie dostępne , faktycznie mogę je zobaczyć i pracować z plikami. gdaldemNarzędziowy gdal można obliczyć nachylenie i hillshade stosując skalę opcję stosunku jednostek pionowej do poziomej. Instrukcja zaleca 111120 m / ° dla czegoś takiego jak płytki SRTM. Na przykład z powłoki OSGeo4W:

$ gdaldem slope -s 111120 -compute_edges N44E007.hgt N44E007_slope.tif

Ta -compute_edgesopcja sprawia, że ​​krawędzie są bardziej płynne, jeśli chcesz zszyć kilka płytek razem. Lub obliczyć kafelki dla dużego regionu. Wadą techniki „skali” jest to, że odległości w kierunkach EW i NS nie są równe, z wyjątkiem równika, więc w przypadku płytek bliżej biegunów mogą występować dziwne błędne interpretacje nachylenia.

Mike T.
źródło
Warto podkreślić ostatni komentarz: jest to słabe rozwiązanie dla punktów znajdujących się w pobliżu równika. To nie jest mała kwestia „dziwnych fałszywych informacji”: wyniki będą rażąco błędne, szczególnie w miejscach bliższych Polakom niż Równik. Dokumentacja dla gdaldemstanów „W przypadku lokalizacji nieopodal równika najlepiej byłoby zmienić projekt siatki za pomocą gdalwarp przed użyciem gdaldem”. Niestety nie zadziała to w przypadku zestawów danych obejmujących kulę ziemską, chyba że podzielisz je na małe części (być może 74 strefy UTM?), Nie wyświetlisz ich, obliczysz zbocza i nie uzyskasz mozaiki wyników.
whuber
7

Krótko mówiąc, nie ma jednego. Z definicji układ współrzędnych oparty na stopniach nie jest rzutowany. W mowie potocznej mówimy, że WGS84 jest projekcją „geograficzną”, ale to nieprawda, dla wygody.

Chyba pamiętam czytanie o oprogramowaniu lub procesie do dokładnej pracy z modelami wysokości w nieprzewidzianej przestrzeni geograficznej, ale nie mogę teraz tego zlokalizować. W każdym razie byłby to eksperymentalny lub zbudowany sam z procesu typu kodu.

Ahhh, znalazł to: Opracowanie globalnego zestawu danych nachylenia do oszacowania występowania osuwisk spowodowanych trzęsieniami ziemi (USGS). Strona 4 dobrze opisuje problem

... długość jednego stopnia zmienia się w zależności od jego szerokości geograficznej. Na równiku blok o jeden stopień na jeden stopień jest rozsądnie kwadratowy po przeliczeniu na jednostki metrów (111 321 metrów w kierunku x o 110 567 metrów w kierunku y ... ale bliżej biegunów odległości w Kierunek x zmniejsza się w zależności od cosinusa szerokości geograficznej, ze względu na zbieżność południków. Większość pakietów GIS, w tym ArcGIS, działa tylko na kwadratowych pikselach, a więc używając współczynnika do dostosowania wymiarów x, y lub z, aby wspólna jednostka nie jest możliwa.

W dalszej części artykułu opisano konkretne obliczenia i narzędzia programowe ( gdal , python , numpy ), których używali do obejścia tego podstawowego problemu. Artykuł nie zawiera kodu, ale jeśli zostanie o to poproszony, może go udostępnić. W każdym razie prawdopodobnie po prostu zapytam, gdzie są wyniki, ponieważ jako USGS prawdopodobnie jest już gdzieś online. :)

matowe wilkie
źródło
1
Sugestia tego artykułu, że azymutalna projekcja w jednakowej odległości może być wykorzystana do obliczania zboczy, jest błędna i błędna. Rzeczywiście da prawidłowe nachylenie w pobliżu początku projekcji, ale one również będą stopniowo coraz mniej dokładne wraz ze wzrostem odległości do początku.
whuber
dzięki za wskazanie tego. Czytelnicy, koniecznie przeczytajcie również gis.stackexchange.com/a/40464/108 , dla zachowania równowagi
Matt Wilkie,
2

Globalne parametry DEM (gdzie większość formuł opiera się na założeniu przestrzeni euklidesowej) można skutecznie wyprowadzić przy użyciu systemu EQUI7 GRID (Bauer-Marschallinger i in. 2014). EQUI7 GRID dzieli świat na 7 obszarów lądowych, wszystkie rzutowane w jednakowo odległym systemie projekcyjnym z minimalną utratą precyzji. Zobacz przykład globalnego DEM przy rozdzielczości 250 m w EQUI7 GRID. Tutaj możesz znaleźć przykładowy kod, który pokazuje, jak uzyskać globalne parametry DEM za pomocą SAGA GIS. Po zakończeniu wyprowadzania parametrów DEM w systemie EQUI7 GRID można przekształcić wszystkie mapy we longlatwspółrzędne WGS84, a następnie utworzyć globalną mozaikę za pomocą GDAL.

Tom Hengl
źródło
Czy możesz wyjaśnić, jak to odpowiada na pytanie? Jeśli proponujesz stosowanie równoodległych rzutów do obliczeń nachylenia, pamiętaj, że jest to złe rozwiązanie z powodu dużych względnych zniekształceń metrycznych, które występują, gdy ktoś odchodzi od środka projekcji. Chociaż skupienie siedmiu takich prognoz na masach lądowych pomaga złagodzić ten problem, wciąż nie jest to najlepszy wybór.
whuber
Artykuł Bauer-Marschallinger i in. (2014) wyjaśnia, dlaczego te prognozy wybrano do reprezentowania globalnych mas lądowych (zakłada się, że mają minimalną utratę precyzji). Zgadzam się, że każda projekcja 2D ostatecznie doprowadzi do deformacji, ale o ile mi wiadomo, EQUI7 jest dobrym kompromisem między utratą precyzji a towarem (algebra 2D). Powiedziawszy to, sześciokąty są ponownie używane do reprezentowania globalnych powierzchni ziemi (chociaż analiza DEM z sześciokątami 3D jest nadal uciążliwa).
Tom Hengl
Dziękuję za referencje. Jej streszczenie sugeruje, że rozwiązuje on zupełnie inny problem, polegający na „minimalizowaniu nadpróbkowania danych lokalnych pojawiających się podczas projekcji ogólnych zdjęć satelitarnych na regularną siatkę rastrową”. Nie oznacza to, że będzie dobrze działać w innych celach, takich jak szacowanie nachyleń.
whuber
Oczywiście EQUI7 nie rozwiązuje absolutnie problemu dokładnego oszacowania lokalnych nachyleń, ale jest to prawdopodobnie bardziej eleganckie rozwiązanie niż użycie wyżej sugerowanej projekcji Mercatora. Ostatecznie, jeśli chcemy oszacować nachylenie z idealną precyzją, wówczas prawdopodobnie jedynymi opcjami są (1) zastosowanie lokalnych (jednakowo odległych) rzutów na płytki o mniejszych rozmiarach (np. 100 na 100 km) z 10-20% zachodzeniem, jak wspomniano również w Verdin i in. (2007) lub (2), aby użyć siatki sześciokątnej ( pakiet dggridR ).
Tom Hengl
Problemem nie jest precyzja - polega ona na tworzeniu systematycznie tendencyjnych nachyleń i aspektów. Ponieważ jednakowo odległe rzuty różnicowo zniekształcają kierunki ortogonalne względem geodezji pochodzącej z ich centrów, aspekty zawsze będą błędne (chociaż dość dokładne w pobliżu centrów, w których wszystkie zniekształcenia są niskie), a błędy na zboczach będą rosły gwałtownie wraz z odległością. Wykorzystanie wielu lokalnych projekcji oczywiście zadziała, ale jest przeciwieństwem elegancji, którą cenisz.
whuber
-2

Nachylenie to wzrost / bieg. Oblicz wzrost i uruchom obliczenia, a otrzymasz odpowiedź. Obliczenie odległości między współrzędnymi geograficznymi jest proste. Spowoduje to zmniejszenie błędu ponownego próbkowania w porównaniu do konwersji na UTM itp.

THK
źródło