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

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.

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.

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)
}

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

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.

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.

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?

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.