Uzyskaj wartości rastrowe z nakładki wielokątów w rozwiązaniach GIS opensource

16

Mam dwie warstwy. Warstwa w kształcie wielokąta z wieloma płytkami i warstwa rastrowa zawierająca pokrycie terenu CORINE 2006 z wieloma kategoriami w mapie kolorów. Chcę uzyskać dla każdego wielokąta w kształcie kapsuły sumę każdej kategorii pokrycia terenu warstwy rastrowej.

Na przykład istnieje wielokąt o identyfikatorze „2” i chcę przypisać takie atrybuty dla tego wielokąta (w procentach lub metrach kwadratowych):

  • Grunty orne: 15%
  • Las: 11%
  • Ulice: 2% (... i tak jeden)

Próbowałem to zrobić na trawie, qgis (bez funkcji), saga (po prostu sumuje każdy do całkowitej wartości) r (całkowita suma), ale nadal nie znalazłem rozwiązania. Większość wtyczek (statystyki strefowe w qgis) obsługuje tylko warstwy rastrowe 0-1. v.rast.stats też nie pomogło. Jestem otwarty na każde dobre i inteligentne rozwiązanie! Może nawet zastosowałem niewłaściwe podejście lub popełniłem błędy.

W Arcgis to zadanie jest dość łatwe, jeśli dobrze pamiętam, ale wciąż brakuje mi dobrego rozwiązania dla twojego codziennego użytkownika linuksa.

Korzystam z systemu Linux Debian i dlatego mogę używać programów tylko dla tego systemu operacyjnego.


EDYCJA: Ponieważ to pytanie wciąż ma tak wiele wyświetleń i odwiedzających: Napisałem wtyczkę QGIS, która jest również w stanie obliczyć pokrycie warstwy rastrowej. Nie zakodowałem jeszcze nakładki wielokąta, ale na pewno jest to zaplanowane. Znajdź wtyczkę tutaj i zainstaluj najpierw bibliotekę Scipy.

Kulik
źródło
Można to zdecydowanie zrobić w R, to tylko kwestia ustalenia, które funkcje. Musisz nałożyć raster na każdy wielokąt, a następnie użyć table (), aby uzyskać podsumowanie pikseli „wycinanych przez pliki cookie”. Przydatne mogą być pakiety raster, rgdal i rgeos. Przeczytaj „Widok przestrzenny zadania R” (google go znajdzie)
Spacedman,
jasne, ale jak mogę uzyskać takie streszczenie. Możesz z łatwością nakładać warstwę wielokąta warstwą rastrową za pomocą! Is.na (nakładka (Poly, Raster)), ale za pomocą poleceń takich jak ekstrakt mogę obliczyć całkowity obszar w pikselu wycinanym ciasteczkami, a nie różne kategorie mapy kolorów . Nie znałem rgeo, ale przejrzałem API i nie znalazłem żadnej funkcji, aby to zrobić.
Curlew,
Sprawdź r.univar w GRASS, jak patrz grasswiki.osgeo.org/wiki/Zonal_statistics
markusN
Cześć! Dziękujemy za zrobienie wtyczki QGIS! Chciałem tylko wspomnieć, że wtyczka ulega awarii (> 13000 wielokątów). Byłoby wspaniale, gdyby podzielił zadanie, aby nie upaść. I byłoby wspaniale mieć opcję dodania wszystkich klas jednocześnie (np. Aby tabela atrybutów otrzymała 2 nowe pola LandcoverID i Landcover%, gdzie oba zawierają listę CSV z wartościami) :)
Mfbaer
@Joran: Jeśli uważasz, że to błąd, zgłoś raport o błędzie zamiast pisać go w komentarzu ( github.com/Martin-Jung/LecoS/issues ). Ponadto 1) serializowanie lub przetwarzanie wsadowe zadań nie jest zadaniem wtyczek. Następnie uruchom go na mniejszych podzbiorach. 2) Jasne. Można by dodać wiele wspaniałych rzeczy. Kod jest open source, możesz go kodować :)
Curlew,

Odpowiedzi:

13

Użyj „wyciągu”, aby nałożyć elementy wielokąta ze SpatialPolygonsDataFrame (który można odczytać z pliku kształtu za pomocą maptools: readShapeSpatial) na rastrze, a następnie użyj „tabeli” do podsumowania. Przykład:

> c=raster("cumbria.tif") # this is my CORINE land use raster
> summary(spd)
Object of class SpatialPolygonsDataFrame
[...]
> nrow(spd)  # how many polygons:
[1] 2
> ovR = extract(c,spd)
> str(ovR)
List of 2
 $ : num [1:542] 26 26 26 26 26 26 26 26 26 26 ...
 $ : num [1:958] 18 18 18 18 18 18 18 18 18 18 ...

Więc mój pierwszy wielokąt obejmuje 542 piksele, a mój drugi obejmuje 958. Mogę podsumować każdy z nich:

> lapply(ovR,table)
[[1]]

 26  27 
287 255 

[[2]]

  2  11  18 
 67  99 792 

Tak więc mój pierwszy wielokąt ma 287 pikseli z klasy 26 i 255 pikseli z klasy 27. Łatwo jest zsumować, podzielić i pomnożyć przez 100, aby uzyskać wartości procentowe.

Spacedman
źródło
Świetnie, wielkie dzięki za wysiłek. Spróbuję i
zdam
6

Chciałem się zgłosić i oto jestem. Rozwiązanie Spacedman działało świetnie i mogłem wyeksportować wszystkie informacje dla każdego wielokąta w moim kształcie. Na wypadek, gdyby ktoś miał ten sam problem, oto jak poprzedziłem:

...
tab <- apply(ovR,table)
# Calculate percentage of landcover types for each polygon-field.
# landcover is a datastream with the names of every polygon
for(i in 1:length(tab)){
 s <- sum(tab[[i]])
 mat <- as.matrix(tab[[i]])
 landcover[i,paste("X",row.names(mat),sep="")] <- as.numeric(tab[[i]]/s)
}
Kulik
źródło
3

jeśli dobrze rozumiem, czego chcesz, i zakładając, że masz warstwę wektorową „mypolygonlayer” i warstwę rastrową „corina” już w bazie danych GRASS GIS:

Najpierw przekonwertowałem wektor na raster. Kot zapewni, że będziesz mieć jeden unikalny identyfikator na wielokąt. Jeśli masz kolumnę z unikalnym identyfikatorem numerycznym, możesz zamiast niej użyć tej kolumny. Kolumna etykiet jest opcjonalna:

v.to.rast input = warstwa mypolygonlayer = 1 wynik = użycie mypolygons = cat labelcolumn = NameMappingUnit

Następnie uruchom r.stats, aby uzyskać statystyki:

r.stats -a -l input = mypolygons, corina separator =; wyjście = / home / paulo / corinastats.csv

Ostatnim krokiem jest otwarcie pliku corinastats.csv np. W LibreOffice i utworzenie tabeli przestawnej lub użycie R do utworzenia tabeli krzyżowej

Ecodiv
źródło
3

Wiem, że ten post jest dość stary, ale ostatnio nie mogę przeprowadzić tego samego rodzaju analizy, ale pobieranie programów takich jak R jest nieco kłopotliwe na moim komputerze i wymaga zatwierdzenia. Po wielu godzinach poszukiwań metody, którą mogłem zastosować tylko z QGis i Excelem, odkryłem, że ta metoda działa najlepiej dla mnie i chciałam ją zaoferować osobom w takiej samej sytuacji.

  1. Przytnij wielokąt do warstwy rastrowej (Raster → ekstrakcja → clipper: plik wejściowy = warstwa rastrowa, wybierz nazwę wyjściową i lokalizację, kliknij warstwę maski, wybierz wielokąt → ok)

  2. Poligonizuj warstwę strzyżenia (Raster → Konwersja → poligonizacja: plik wejściowy = twoja warstwa klipu, zapisz dane wyjściowe → ok)

  3. Obliczanie liczby pikseli (kliknij właśnie utworzony plik kształtu → otwórz kalkulator pola: zaznacz „utwórz nowe pole” i dodaj nazwę pola, Funkcja = geometria → obszar → ok). Powinieneś teraz mieć nową kolumnę w tabeli atrybutów pokazującą liczbę pikseli.

  4. Zapisz warstwę poligoniczną (kliknij warstwę prawym przyciskiem myszy, zapisz jako: format = plik DBF, zapisz jako → ok)

  5. Podsumowując liczbę pikseli dla każdego siedliska (uruchom Excel, otwórz plik, jeśli nie masz tytułów, dodaj teraz jeden dla każdej kolumny, kliknij pustą komórkę, przejdź do karty DANE, skonsoliduj, upewnij się, że jest suma, kliknij czerwona strzałka obok „Przeglądaj…” i wybierz dwie kolumny (w tym tytuły), kliknij „dodaj” i zaznacz oba pola „Górny wiersz” i „lewa kolumna” → ok)

  6. Jeśli, podobnie jak ja, masz wiele wielokątów do analizy i musisz je porównać w tej samej tabeli, ten krok będzie pomocny. W nowym skoroszycie Excela wypisz numery siedlisk w kolumnie A (dla mnie od 1 do 48) i umieść dwie kolumny, które właśnie skonsolidowałeś w kolumnach B i C (siedlisko w B i liczba pikseli w C). W komórce D1 napisz następującą formułę: = IFNA (INDEKS (C: C; MECZ (A2; B: B; 0)); "") i przeciągnij (lub dwukrotnie kliknij prawy dolny róg) do ostatniej wartości (więc jeśli masz 48 siedlisk aż do komórki D48). Liczba pikseli znajduje się teraz w odpowiednich komórkach siedliska i możesz powtórzyć ten proces dla wszystkich swoich wielokątów.

Emilia
źródło
2

Co powiesz na konwersję danych CORINE na zbiór danych wielokątów wektorowych za pomocą QGIS ( Raster> Konwersja> Poligonizacja ), a następnie za pomocą funkcji Union ( Vector> Narzędzia geoprzetwarzania> Zjednoczenie ) w celu połączenia z wielokątami. Wynikowy zestaw danych wektorowych zawierałby obszary każdej klasy CORINE w każdym wielokącie.

Andy Harfoot
źródło
dzięki za tę sugestię. Nie myślałem jeszcze o unii wektorowej. Może spróbuję tego, jeśli przetwarzanie R w jakiś sposób zawiedzie.
Curlew
0

QGIS.

W bagażniku QGIS dostępna jest inna wersja ZonalStats, nazywa się to Zonal Statistics.

Wykonuje to wymaganą funkcję.

Jeśli chodzi o przepływ pracy, nie jestem pewien, ile masz rastrów, czy są to tylko zespoły rastrowe?

Willy
źródło
dzięki za komentarz, ale statystyki strefowe jedzą tylko raster bez kategorii. Korzystam z QGIS Trunk 1.9
Curlew,
0

W przeciwieństwie do większości powyższych odpowiedzi, argumentowałbym, że lepszym rozwiązaniem jest rasteryzacja twoich wielokątów i praca z dwoma zestawami danych rastrowych zamiast z dwoma zestawami danych wielokątów. Jest to znacznie mniej intensywne przetwarzanie i dlatego jest jedynym łatwym do wdrożenia rozwiązaniem do przetwarzania dużych plików rastrowych i dużych plików wielokątów w R.

Po zrasteryzowaniu wielokątów w dokładnie tym samym zakresie i rozdzielczości danych rastrowych, możesz tabelarycznie podsumować statystyki, jak wyjaśniono tutaj , co jest odpowiednie, jeśli twój raster mieści się w pamięci (małe / średnie warstwy rastrowe) lub możesz binaryzować każdą kategorię za pomocą reclassfunkcji i niż obliczyć zonalstatystyki dla każdej klasy. Oto rozwiązanie, które łączy rasteryzację i statystyki strefowe w jedną funkcję i działa dobrze z bardzo dużymi zestawami danych.

joaoal
źródło