Selekcja i modyfikacja danych
Wstęp
Selekcja podzbiorów danych jest jednym z najważniejszych działań programistycznych, a bazy danych są częścią praktycznie każdej aplikacji. W Systemach Informacji Geograficznej selekcja podzbiorów danych jeną z najczęściej podejmowanych czynności. W SIG rozróżniamy dwa typy zapytań: zapytania atrybutowe, niekiedy określane jako nieprzestrzenne, analizujące wartości atrybutów oraz bardziej złożone zapytania przestrzenne, wykorzystujące różne relacje przestrzenne. R poprzez bibliotekę rgeos pozwala wykonywać zarówno zapytania przestrzenne jak i atrybutowe.
Wymagane biblioteki: rgeos, raster, sp
Wymagana wiedza: systemy baz danych, podstawy systemów Inf. , przestrzenne bazy danych
Wczytanie danych wejściowych
Z bazy danych geometry.sqlite interesują nas warstwy: obszar, parcele, droga, trasa, odcinki, point1, point2
library(RSQLite)## Loading required package: DBIconn <- dbConnect(SQLite(),"") #obiekt połączenie
dbListTables(conn) # funkcja wyświetla listę tabel## character(0)dbDisconnect(conn)## [1] TRUEwarstwy <- c("obszar", "parcele", "droga", "trasa", "odcinki", "point1", "point2")
typy <- c("wkbPoints","wkbLineStrings","wkbPolygons")
library(rgdal)## Loading required package: sp## rgdal: version: 1.1-10, (SVN revision 622)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 2.1.0, released 2016/04/25
## Path to GDAL shared files: /usr/share/gdal/2.1
## Loaded PROJ.4 runtime: Rel. 4.9.2, 08 September 2015, [PJ_VERSION: 492]
## Path to PROJ.4 shared files: (autodetected)
## Linking to sp version: 1.2-3#sprawdzenie typów
sapply(warstwy, function(layer,dsn) {ogrInfo(dsn=dsn,layer=layer)$eType},dsn="DATA/geometry.sqlite")## obszar parcele droga trasa odcinki point1 point2
## 3 3 2 2 2 1 1Selekcja danych wektorowych
Selekcja na podstawie atrybutów
Selekcja na podstawie atrybutów jest naturalnym sposobem selekcji, i przypomina konstruowanie zapytań w języku SQL. R pozwala również na selekcję poprzez wskazanie indeksów. W obu przypadkach selekcja wykonywana jest zgodnie z zasadami selekcji wierszy w data.frame. W przypadku Spatial*DataFrame wybierane są zarówno atrybuty ze slotu data jak i odpowiednie elementy geometrii. Należy tu jednak zaznaczyć, że istotą R jest wektoryzacja, co oznacza że proces selekcji nie polega na testowania poszczególnych wierszy czy spełniają kryteria wyszukiwania ale na przygotowaniu wektora logicznego, o długości zbioru danych i wyboru tych elementów zbioru, które spełniają kryterium TRUE.
dns <- "DATA/geometry.sqlite"
punkty <- readOGR(dns,"point2")## OGR data source with driver: SQLite
## Source: "DATA/geometry.sqlite", layer: "point2"
## with 22 features
## It has 3 fieldspoligony <- readOGR(dns,"parcele")## OGR data source with driver: SQLite
## Source: "DATA/geometry.sqlite", layer: "parcele"
## with 9 features
## It has 2 fieldsplot(punkty)
#Selekcja poprzez wskazanie indeksów wierszy
sel <- punkty[c(1,2,3),]
str(sel)## Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
## ..@ data :'data.frame': 3 obs. of 3 variables:
## .. ..$ id : int [1:3] 2 3 2
## .. ..$ value: num [1:3] 1 14.2 11.4
## .. ..$ nazwa: Factor w/ 19 levels "A","A1","A2",..: 1 7 8
## ..@ coords.nrs : num(0)
## ..@ coords : num [1:3, 1:2] 356176 356345 356501 507031 506960 ...
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : NULL
## .. .. ..$ : chr [1:2] "coords.x1" "coords.x2"
## ..@ bbox : num [1:2, 1:2] 356176 506819 356501 507031
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : chr [1:2] "coords.x1" "coords.x2"
## .. .. ..$ : chr [1:2] "min" "max"
## ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
## .. .. ..@ projargs: chr "+proj=tmerc +lat_0=0 +lon_0=19 +k=0.9993 +x_0=500000 +y_0=-5300000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"plot(punkty)+points(sel,col=2,pch=1,lwd=2,cex=2)
## numeric(0)#Zapytanie w stylu SQL
selector <- punkty$id==2 & punkty$value>0 #przygotowanie wektora selekcji
selector## [1] TRUE FALSE TRUE TRUE FALSE FALSE FALSE FALSE TRUE FALSE FALSE
## [12] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSEsel <- punkty[selector,] #zastosowanie wektora selekcji,
sel <- punkty[punkty$id==2 & punkty$value>0,] #równoważne, ale mniej czytelne
plot(punkty)+points(sel,col=2,pch=1,lwd=2,cex=2)
## numeric(0)#zapytanie dla obiektu zawierającego poligony
sel <- poligony[grep("^A",poligony$nazwa),] # grep pozwala na częściowe dopasowanie łańcucha tekstu na podstawie wyrażenia regularnego.
plot(poligony) + plot(sel,col=2,lwd=2,add=TRUE)
## numeric(0)Wektoryzacja powoduje, że wektor selekcji nie musi być związany przetwarzanym zbiorem, ale może być od niego niezależny. Na przykład zbiór losowy.
selector <- sort(sample(1:length(punkty),10)) # 10 punktów losowych
sel <- punkty[selector,]
plot(punkty)+points(sel,col=2,pch=1,lwd=2,cex=2)
## numeric(0)Selekcja na podstawie współrzędnych
Selekcja na podstawie współrzędnych może być w praktyce stosowana to prostych typów danych: SpatialPoints- SpatialGrid- oraz SpatialGrid- DataFrame oraz dla obiektów typu Raster*. Dla typów Polygon i LineString, ze względu na złożoność struktury wewnętrznej takie zapytanie jest dużo trudniejsze i powinno być zastępowane zapytaniami opartymi o relacje przestrzenne.
Wektor
Przykład zapytania, które wybiera punkty o współrzędnych powyżej wsp.x i wsp.y. Podobnie jak w przypadku selekcji atrybutowej wykorzystana jest wektoryzacja. Tworzony jest wektor selekcji a następnie stosowany na wybranej *DataFrame. Zadaniem programisty jest prawidłowe przygotowanie kryterium selekcji, na przykład na podstawie zmiennych.
wsp.x <- 356161
wsp.y <- 506420
selector <- coordinates(punkty)[,1]>wsp.x & coordinates(punkty)[,2]>wsp.y
sel <- punkty[selector,]
plot(punkty)+points(sel,col=2,pch=1,lwd=2,cex=2)
## numeric(0)Kryteria selekcji mogą być wyliczone z formuł matematycznych, na przykład jako obszar dookoła środka zbioru danych \(\pm\) odchylenie standardowe. W takiej sytuacji selektor jest wynikiem kilku operacji logicznych i zdecydowanie lepiej jest przygotować osobny wektor.
cr.x <- coordinates(punkty)[,1]
cr.y <- coordinates(punkty)[,2]
sel.x <- cr.x < (mean(cr.x) + sd(cr.x)) & cr.x > (mean(cr.x) - sd(cr.x))
sel.y <- cr.y < (mean(cr.y) + sd(cr.y)) & cr.y > (mean(cr.y) - sd(cr.y))
selector <- sel.x & sel.y #operacja logiczna
selector## [1] FALSE FALSE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE FALSE
## [12] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSEsel <- punkty[selector,]
plot(punkty)+points(sel,col=2,pch=1,lwd=2,cex=2)
## numeric(0)Selekcja na podstawie realcji przestrzennych
Relacje i operatory przestrzenne pozwalają nam wybrać obiekty z jednej warstwy na podstawie relacji przestrzennej względem obiektu(ów) z innej warstwy. Część relacji ma charakter binarny (spełnia kryterium lub nie) do takich relacji należy na przykład within, overlap, contains itp. Druga część operacji opiera się na operatorach związanych z odległością pomiędzy obiektami lub ich powierzchnią, długością i obwodem. Szczegółowe omówienie wszystkich możliwych relacji na obecnym etapie wykracza poza zakres skryptu, po szczegóły proponuję zapoznać się z dokumentacją rgeos. Zastosowanie relacji omówimy na podstawie kilku przykładów
Wybierz poligony będąc w relacji do innego poligonu
library(rgeos)## rgeos version: 0.3-21, (SVN revision 540)
## GEOS runtime version: 3.5.0-CAPI-1.9.0 r4084
## Linking to sp version: 1.2-3
## Polygon checking: TRUEdns <- "DATA/geometry.sqlite"
pselektor <- readOGR(dns,"obszar") # wybieramy selektor## OGR data source with driver: SQLite
## Source: "DATA/geometry.sqlite", layer: "obszar"
## with 2 features
## It has 1 fieldspoligony <- readOGR(dns,"parcele")## OGR data source with driver: SQLite
## Source: "DATA/geometry.sqlite", layer: "parcele"
## with 9 features
## It has 2 fieldsplot(pselektor,col="yellow")+plot(poligony, add=TRUE)
## numeric(0)str(pselektor)## Formal class 'SpatialPolygonsDataFrame' [package "sp"] with 5 slots
## ..@ data :'data.frame': 2 obs. of 1 variable:
## .. ..$ id: int [1:2] 1 2
## ..@ polygons :List of 2
## .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
## .. .. .. ..@ Polygons :List of 1
## .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
## .. .. .. .. .. .. ..@ labpt : num [1:2] 357020 506688
## .. .. .. .. .. .. ..@ area : num 2155287
## .. .. .. .. .. .. ..@ hole : logi FALSE
## .. .. .. .. .. .. ..@ ringDir: int 1
## .. .. .. .. .. .. ..@ coords : num [1:12, 1:2] 357077 356564 356374 356279 356500 ...
## .. .. .. ..@ plotOrder: int 1
## .. .. .. ..@ labpt : num [1:2] 357020 506688
## .. .. .. ..@ ID : chr "1"
## .. .. .. ..@ area : num 2155287
## .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
## .. .. .. ..@ Polygons :List of 1
## .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
## .. .. .. .. .. .. ..@ labpt : num [1:2] 359010 506369
## .. .. .. .. .. .. ..@ area : num 616287
## .. .. .. .. .. .. ..@ hole : logi FALSE
## .. .. .. .. .. .. ..@ ringDir: int 1
## .. .. .. .. .. .. ..@ coords : num [1:11, 1:2] 358899 358682 358702 358757 358954 ...
## .. .. .. ..@ plotOrder: int 1
## .. .. .. ..@ labpt : num [1:2] 359010 506369
## .. .. .. ..@ ID : chr "2"
## .. .. .. ..@ area : num 616287
## ..@ plotOrder : int [1:2] 1 2
## ..@ bbox : num [1:2, 1:2] 356279 505612 359340 507692
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : chr [1:2] "x" "y"
## .. .. ..$ : chr [1:2] "min" "max"
## ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
## .. .. ..@ projargs: chr "+proj=tmerc +lat_0=0 +lon_0=19 +k=0.9993 +x_0=500000 +y_0=-5300000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"#są dwa poligony, wybieramy poligon nr1
pol_selector <- pselektor[pselektor$id==1,]
#Wybierz poligony zawarte w całości w innym poligonie
selector <- gWithin(poligony,pol_selector,byid = TRUE)
str(selector) #zwraca macierz, w tym wypadku macierz ma tylko 1 wiersz, dlatego wybieramy [1,]## logi [1, 1:9] FALSE TRUE FALSE FALSE FALSE TRUE ...
## - attr(*, "dimnames")=List of 2
## ..$ : chr "1"
## ..$ : chr [1:9] "1" "2" "3" "4" ...plot(poligony) + plot(poligony[selector[1,],],lwd=2,col="yellow",add=TRUE)## numeric(0)selector <- gContains(pol_selector,poligony,byid = TRUE) #contains to odwrotność within
str(selector) # tym razem macierz jest odwrotna wybieramy [,1]## logi [1:9, 1] FALSE TRUE FALSE FALSE FALSE TRUE ...
## - attr(*, "dimnames")=List of 2
## ..$ : chr [1:9] "1" "2" "3" "4" ...
## ..$ : chr "1"plot(poligony) + plot(poligony[selector[,1],],lwd=2,col="yellow",add=TRUE)
## numeric(0)#Wybierz poligony zawierające się lub przecinające się z innym poligonem
selector <- gIntersects(pol_selector,poligony,byid = TRUE) #contains to odwrotność within
str(selector) # tym razem macierz jest odwrotna wybieramy [,1]## logi [1:9, 1] TRUE TRUE TRUE FALSE TRUE TRUE ...
## - attr(*, "dimnames")=List of 2
## ..$ : chr [1:9] "1" "2" "3" "4" ...
## ..$ : chr "1"plot(poligony) + plot(poligony[selector[,1],],lwd=2,col="yellow",add=TRUE)
## numeric(0)#Wybierz poligony tylko przecinające się z innym poligonem
selector <- gOverlaps(pol_selector,poligony,byid = TRUE)
str(selector) # tym razem macierz jest odwrotna wybieramy [,1]## logi [1:9, 1] TRUE FALSE TRUE FALSE TRUE FALSE ...
## - attr(*, "dimnames")=List of 2
## ..$ : chr [1:9] "1" "2" "3" "4" ...
## ..$ : chr "1"plot(poligony) + plot(poligony[selector[,1],],lwd=2,col="yellow",add=TRUE)
## numeric(0)Wybierz poligony poprzez zasięg przestrzenny
Aby uniknąć skomplikowanych operacji na zagnieżdżonych listach zastosujemy kilka kroków, które przekształcą 4 wartości wyznaczające oczekiwany (lub otrzymany w wyniku obliczeń) zasięg przestrzenny na SpatialPolygon, a następnie tak otrzymany SpatialPolygon użyjemy jako narzędzie selekcji. W tym celu dokonamy konwersji na obiekt extent pakietu raster, a następnie ten obiekt zamienimy funkcją as()
library(raster)
plot(poligony)
xmin <- 356327
xmax <- 357841
ymin <- 506261
ymax <- 507469
ext <- extent(xmin,xmax,ymin,ymax) # budujemy extent
poligon_ext <- as(ext,"SpatialPolygons") # i konwertujemy go na spatialPolygons (nie dataFrame!)
crs(poligon_ext) <- crs(poligony) # musimy ujednolicić georeferencję
selector <- gContains(poligon_ext,poligony,byid = TRUE) # tylko całe a więc gContains
str(selector) ## logi [1:9, 1] FALSE TRUE FALSE FALSE FALSE FALSE ...
## - attr(*, "dimnames")=List of 2
## ..$ : chr [1:9] "1" "2" "3" "4" ...
## ..$ : chr "1"plot(poligony) + plot(poligony[selector[,1],],lwd=2,col="yellow",add=TRUE)
## numeric(0)Może się okazać, ze zasięg przestrzenny otrzymamy jako listę współrzędnych punktów. Mogą to być współrzędne wprowadzone przez użytkownika skryptu jako parametr, mogą też być to współrzędne będące efektem dowolnych operacji geo-przestrzennych. W takiej sytuacji, aby uniknąć skomplikowanych operacji budowania spatialPolygon, możemy skorzystać z polecenia spPolygons().
#współrzędne trzech punktów:
p1 <- c(356512,507741)
p2 <- c(357515,506553)
p3 <- c(356288,506660)
p <- spPolygons(rbind(p1,p2,p3,p1)) # poligon wymaga zamknięcia pierwszym punktem
crs(p) <- crs(poligony)
selector <- gOverlaps(p,poligony,byid = TRUE) # wybieramy te które są całe lub się przecinają
str(selector) ## logi [1:9, 1] TRUE TRUE FALSE FALSE TRUE TRUE ...
## - attr(*, "dimnames")=List of 2
## ..$ : chr [1:9] "1" "2" "3" "4" ...
## ..$ : chr "1"plot(poligony)+plot(p,border=2,lwd=2,add=TRUE)+ plot(poligony[selector[,1],],lwd=2,col="yellow",add=TRUE)
## numeric(0)Selekcja obiektów na podstawie operatora odległości
W przypadku obiektów punktowych proste operacje binarne rzadko zdają egzamin. Zdecydowanie prościej dokonywać analizy poprzez operator odległości, który pozwala dokonać selekcji, jeżeli obiekty znajdują się w pewnej odległości od innego obiektu.
library(rgeos)
library(rgdal)
dns <- "DATA/geometry.sqlite"
linia <- readOGR(dns,"trasa") # wybieramy selektor## OGR data source with driver: SQLite
## Source: "DATA/geometry.sqlite", layer: "trasa"
## with 1 features
## It has 1 fieldspunkty <- readOGR(dns,"point2")## OGR data source with driver: SQLite
## Source: "DATA/geometry.sqlite", layer: "point2"
## with 22 features
## It has 3 fieldsstr(punkty) # obiekt punkty ma trzy atrybuty## Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
## ..@ data :'data.frame': 22 obs. of 3 variables:
## .. ..$ id : int [1:22] 2 3 2 2 3 4 2 2 2 3 ...
## .. ..$ value: num [1:22] 1 14.2 11.4 7 14 ...
## .. ..$ nazwa: Factor w/ 19 levels "A","A1","A2",..: 1 7 8 9 11 14 15 16 17 18 ...
## ..@ coords.nrs : num(0)
## ..@ coords : num [1:22, 1:2] 356176 356345 356501 356405 356483 ...
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : NULL
## .. .. ..$ : chr [1:2] "coords.x1" "coords.x2"
## ..@ bbox : num [1:2, 1:2] 355874 505639 358152 507143
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : chr [1:2] "coords.x1" "coords.x2"
## .. .. ..$ : chr [1:2] "min" "max"
## ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
## .. .. ..@ projargs: chr "+proj=tmerc +lat_0=0 +lon_0=19 +k=0.9993 +x_0=500000 +y_0=-5300000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"dst <- gDistance(linia,punkty,byid = TRUE) #obiekt (macierz z 1 kolumną) dst posłuży do selekcji
selector <- dst[,1]<300 & punkty$value<10 #dwa atrybuty, jeden przestrzenny (dst) drugi nieprzestrzenny (atrybut value)
sel <- punkty[selector,] #
plot(punkty)+plot(linia,add=TRUE)+points(sel,col=2,pch=1,lwd=2,cex=2)
## numeric(0)Selekcja obiektów na podstawie atrybutów z innej warstwy
Narzędzia selekcji wbudowane w biblioteki i środowisko R (i właściwe dla wszystkich zaawansowanych środowisk GIS) poprzez operacje przestrzenne pozwalają wybierać obiekty z jednej warstwy, na podstawie atrybutów innej warstwy. Rozważmy przypadek, przedstawiony na poniższej rycinie. Zależy nam na wyborze poligonów, w których znajdują się punkty z wartością value większą niż 15 (np. parcele zawierające drzewa starsze niż 15 lat)
dns <- "DATA/geometry.sqlite"
poligony <- readOGR(dns,"parcele") # wybieramy selektor## OGR data source with driver: SQLite
## Source: "DATA/geometry.sqlite", layer: "parcele"
## with 9 features
## It has 2 fieldspunkty <- readOGR(dns,"point2")## OGR data source with driver: SQLite
## Source: "DATA/geometry.sqlite", layer: "point2"
## with 22 features
## It has 3 fieldsplot(poligony)+plot(punkty,add=TRUE)+plot(punkty[punkty$value>15,],pch=1,col="green",cex=2,add=TRUE)
## numeric(0)Aby wybrać te punkty, należy przeprowadzić dwie operacje: 1. operację selekcji atrybutowej (nieprzestrzennej), która jest prostsza i mniej wymagająca obliczeniowo 2. operację selekcji przestrzennej na podstawie wyników wcześniejszej operacji.
x <- 15 #zmienna, przekazywana jako parametr
punkty_sel <- punkty[punkty$value>x,] #prosta selekcja
selector <- gContains(poligony,punkty_sel,byid=TRUE)
#lub w jednym zapytaniu. mniej czytelne
selector <- gContains(poligony,subset(punkty,value>x),byid=TRUE) # do selekcji możemy zastosować funkcję subset
selector # jest to macierz ilość poligonów x ilość punktów## 1 2 3 4 5 6 7 8 9
## 9 FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE
## 11 FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## 19 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## 20 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## 21 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSEPonieważ wynikiem operacji gContains, z parametrem byid = TRUE jest macierz pokazująca relację każdego poligonu do każdego punktu w selekcji a nas interesują poligony, które zawierają przynajmniej jeden element należy “spłaszczyć” macierz do pojedynczego wektora (by columns). Służy do tego funkcja any(), która zwraca prawdę jeżeli przynajmniej jeden warunek jest prawdą. Jej przeciwieństwem jest funkcja all(), która wymaga aby wszystkie obiekty były prawdą.
apply(selector,2,"any") # stosujemy funkcję logiczną any, która## 1 2 3 4 5 6 7 8 9
## FALSE TRUE FALSE FALSE FALSE TRUE FALSE FALSE FALSEselector <- apply(selector,2,"any")
sel <- poligony[selector,]
plot(poligony)+plot(sel,col="yellow",lwd=2,add=TRUE)+plot(punkty[punkty$value>x,],add=TRUE)
## numeric(0)Ponieważ taka selekcja wymaga kilku operacji. można przekształcić je do postaci funkcji:
select_poly_by_val <- function(poligony,punkty,kolumna,x) {
selector <- gContains(poligony,punkty[punkty@data[,kolumna]>x,],byid=TRUE) # kolumna przekazywana jako zmienna, nie możemy zastosować subset
selector <- apply(selector,2,"any")
sel <- poligony[selector,]
sel
}
sel <- select_poly_by_val(poligony,punkty,"value",15)
plot(poligony)+plot(sel,col="yellow",lwd=2,add=TRUE)+plot(punkty[punkty$value>x,],add=TRUE)
## numeric(0)Modyfikacja warstw wektorowych przez inną warstwę
Oprócz operatorów logicznych można stosować funkcje, które na podstawie operacji logicznych na zbiorach przekształcają dane wektorowe. Są to: różnica, suma (unia), przecięcie (intersekcja) oraz różnica symetryczna. Operacje te można przeprowadzić zarówno w bibliotece raster, gdzie można dokonywać przekształceń obiektów rastrowych, wektorowych oraz jednych i drugich oraz rgeos, gdzie operacje ograniczają się jedynie do obiektów wektorowych. Jednakże operacje z biblioteki raster są uproszczone i wielu efektów z biblioteki rgeos, zwłaszcza w sytuacji, gdy mieszamy geometrię nie jest możliwe.
library(raster)
library(rgeos)
dns <- "DATA/geometry.sqlite"
poligon1 <- readOGR(dns,"obszar") # wybieramy selektor## OGR data source with driver: SQLite
## Source: "DATA/geometry.sqlite", layer: "obszar"
## with 2 features
## It has 1 fieldspoligon2 <- readOGR(dns,"parcele")## OGR data source with driver: SQLite
## Source: "DATA/geometry.sqlite", layer: "parcele"
## with 9 features
## It has 2 fieldslinie <- readOGR(dns,"droga")## OGR data source with driver: SQLite
## Source: "DATA/geometry.sqlite", layer: "droga"
## with 2 features
## It has 2 fieldspunkty <- readOGR(dns,"point2")## OGR data source with driver: SQLite
## Source: "DATA/geometry.sqlite", layer: "point2"
## with 22 features
## It has 3 fieldsinter_p1_p2 <- gIntersection(poligon1,poligon2,byid=FALSE) # intersection z biblioteki rgeos
str(inter_p1_p2) # zwraca spatialPolygons## Formal class 'SpatialPolygons' [package "sp"] with 4 slots
## ..@ polygons :List of 1
## .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
## .. .. .. ..@ Polygons :List of 7
## .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
## .. .. .. .. .. .. ..@ labpt : num [1:2] 357146 505951
## .. .. .. .. .. .. ..@ area : num 326082
## .. .. .. .. .. .. ..@ hole : logi FALSE
## .. .. .. .. .. .. ..@ ringDir: int 1
## .. .. .. .. .. .. ..@ coords : num [1:6, 1:2] 356951 356573 357117 357586 357518 ...
## .. .. .. .. .. .. .. ..- attr(*, "dimnames")=List of 2
## .. .. .. .. .. .. .. .. ..$ : NULL
## .. .. .. .. .. .. .. .. ..$ : chr [1:2] "x" "y"
## .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
## .. .. .. .. .. .. ..@ labpt : num [1:2] 356537 506466
## .. .. .. .. .. .. ..@ area : num 94710
## .. .. .. .. .. .. ..@ hole : logi FALSE
## .. .. .. .. .. .. ..@ ringDir: int 1
## .. .. .. .. .. .. ..@ coords : num [1:7, 1:2] 356485 356374 356362 356563 356729 ...
## .. .. .. .. .. .. .. ..- attr(*, "dimnames")=List of 2
## .. .. .. .. .. .. .. .. ..$ : NULL
## .. .. .. .. .. .. .. .. ..$ : chr [1:2] "x" "y"
## .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
## .. .. .. .. .. .. ..@ labpt : num [1:2] 356423 506962
## .. .. .. .. .. .. ..@ area : num 95590
## .. .. .. .. .. .. ..@ hole : logi FALSE
## .. .. .. .. .. .. ..@ ringDir: int 1
## .. .. .. .. .. .. ..@ coords : num [1:7, 1:2] 356340 356279 356310 356456 356605 ...
## .. .. .. .. .. .. .. ..- attr(*, "dimnames")=List of 2
## .. .. .. .. .. .. .. .. ..$ : NULL
## .. .. .. .. .. .. .. .. ..$ : chr [1:2] "x" "y"
## .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
## .. .. .. .. .. .. ..@ labpt : num [1:2] 357436 507141
## .. .. .. .. .. .. ..@ area : num 154212
## .. .. .. .. .. .. ..@ hole : logi FALSE
## .. .. .. .. .. .. ..@ ringDir: int 1
## .. .. .. .. .. .. ..@ coords : num [1:6, 1:2] 357529 357630 357633 357337 357209 ...
## .. .. .. .. .. .. .. ..- attr(*, "dimnames")=List of 2
## .. .. .. .. .. .. .. .. ..$ : NULL
## .. .. .. .. .. .. .. .. ..$ : chr [1:2] "x" "y"
## .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
## .. .. .. .. .. .. ..@ labpt : num [1:2] 357513 506555
## .. .. .. .. .. .. ..@ area : num 99957
## .. .. .. .. .. .. ..@ hole : logi FALSE
## .. .. .. .. .. .. ..@ ringDir: int 1
## .. .. .. .. .. .. ..@ coords : num [1:6, 1:2] 357635 357640 357619 357406 357370 ...
## .. .. .. .. .. .. .. ..- attr(*, "dimnames")=List of 2
## .. .. .. .. .. .. .. .. ..$ : NULL
## .. .. .. .. .. .. .. .. ..$ : chr [1:2] "x" "y"
## .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
## .. .. .. .. .. .. ..@ labpt : num [1:2] 356893 507070
## .. .. .. .. .. .. ..@ area : num 183293
## .. .. .. .. .. .. ..@ hole : logi FALSE
## .. .. .. .. .. .. ..@ ringDir: int 1
## .. .. .. .. .. .. ..@ coords : num [1:5, 1:2] 356599 357018 357165 356786 356599 ...
## .. .. .. .. .. .. .. ..- attr(*, "dimnames")=List of 2
## .. .. .. .. .. .. .. .. ..$ : NULL
## .. .. .. .. .. .. .. .. ..$ : chr [1:2] "x" "y"
## .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
## .. .. .. .. .. .. ..@ labpt : num [1:2] 357005 506487
## .. .. .. .. .. .. ..@ area : num 184630
## .. .. .. .. .. .. ..@ hole : logi FALSE
## .. .. .. .. .. .. ..@ ringDir: int 1
## .. .. .. .. .. .. ..@ coords : num [1:7, 1:2] 356715 357153 357253 356768 356791 ...
## .. .. .. .. .. .. .. ..- attr(*, "dimnames")=List of 2
## .. .. .. .. .. .. .. .. ..$ : NULL
## .. .. .. .. .. .. .. .. ..$ : chr [1:2] "x" "y"
## .. .. .. ..@ plotOrder: int [1:7] 1 7 6 4 5 3 2
## .. .. .. ..@ labpt : num [1:2] 357146 505951
## .. .. .. ..@ ID : chr "1"
## .. .. .. ..@ area : num 1138475
## ..@ plotOrder : int 1
## ..@ bbox : num [1:2, 1:2] 356279 505712 357640 507383
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : chr [1:2] "x" "y"
## .. .. ..$ : chr [1:2] "min" "max"
## ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
## .. .. ..@ projargs: chr "+proj=tmerc +lat_0=0 +lon_0=19 +k=0.9993 +x_0=500000 +y_0=-5300000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"plot(inter_p1_p2)
inter_p1_p2 <- intersect(poligon1,poligon2) # intersect z biblioteki raster
str(inter_p1_p2) #zwraca spatialPolygonsDataFrame## Formal class 'SpatialPolygonsDataFrame' [package "sp"] with 5 slots
## ..@ data :'data.frame': 7 obs. of 3 variables:
## .. ..$ id : int [1:7] 1 1 1 1 1 1 1
## .. ..$ ID : int [1:7] 1 2 3 5 6 7 9
## .. ..$ nazwa: Factor w/ 9 levels "Adam","Agata",..: 5 7 2 4 9 1 3
## ..@ polygons :List of 7
## .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
## .. .. .. ..@ Polygons :List of 1
## .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
## .. .. .. .. .. .. ..@ labpt : num [1:2] 356423 506962
## .. .. .. .. .. .. ..@ area : num 95590
## .. .. .. .. .. .. ..@ hole : logi FALSE
## .. .. .. .. .. .. ..@ ringDir: int 1
## .. .. .. .. .. .. ..@ coords : num [1:7, 1:2] 356340 356279 356310 356456 356605 ...
## .. .. .. .. .. .. .. ..- attr(*, "dimnames")=List of 2
## .. .. .. .. .. .. .. .. ..$ : NULL
## .. .. .. .. .. .. .. .. ..$ : chr [1:2] "x" "y"
## .. .. .. ..@ plotOrder: int 1
## .. .. .. ..@ labpt : num [1:2] 356423 506962
## .. .. .. ..@ ID : chr "1"
## .. .. .. ..@ area : num 95590
## .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
## .. .. .. ..@ Polygons :List of 1
## .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
## .. .. .. .. .. .. ..@ labpt : num [1:2] 356893 507070
## .. .. .. .. .. .. ..@ area : num 183293
## .. .. .. .. .. .. ..@ hole : logi FALSE
## .. .. .. .. .. .. ..@ ringDir: int 1
## .. .. .. .. .. .. ..@ coords : num [1:5, 1:2] 356599 357018 357165 356786 356599 ...
## .. .. .. .. .. .. .. ..- attr(*, "dimnames")=List of 2
## .. .. .. .. .. .. .. .. ..$ : NULL
## .. .. .. .. .. .. .. .. ..$ : chr [1:2] "x" "y"
## .. .. .. ..@ plotOrder: int 1
## .. .. .. ..@ labpt : num [1:2] 356893 507070
## .. .. .. ..@ ID : chr "2"
## .. .. .. ..@ area : num 183293
## .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
## .. .. .. ..@ Polygons :List of 1
## .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
## .. .. .. .. .. .. ..@ labpt : num [1:2] 357436 507141
## .. .. .. .. .. .. ..@ area : num 154212
## .. .. .. .. .. .. ..@ hole : logi FALSE
## .. .. .. .. .. .. ..@ ringDir: int 1
## .. .. .. .. .. .. ..@ coords : num [1:6, 1:2] 357529 357630 357633 357337 357209 ...
## .. .. .. .. .. .. .. ..- attr(*, "dimnames")=List of 2
## .. .. .. .. .. .. .. .. ..$ : NULL
## .. .. .. .. .. .. .. .. ..$ : chr [1:2] "x" "y"
## .. .. .. ..@ plotOrder: int 1
## .. .. .. ..@ labpt : num [1:2] 357436 507141
## .. .. .. ..@ ID : chr "3"
## .. .. .. ..@ area : num 154212
## .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
## .. .. .. ..@ Polygons :List of 1
## .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
## .. .. .. .. .. .. ..@ labpt : num [1:2] 356537 506466
## .. .. .. .. .. .. ..@ area : num 94710
## .. .. .. .. .. .. ..@ hole : logi FALSE
## .. .. .. .. .. .. ..@ ringDir: int 1
## .. .. .. .. .. .. ..@ coords : num [1:7, 1:2] 356485 356374 356362 356563 356729 ...
## .. .. .. .. .. .. .. ..- attr(*, "dimnames")=List of 2
## .. .. .. .. .. .. .. .. ..$ : NULL
## .. .. .. .. .. .. .. .. ..$ : chr [1:2] "x" "y"
## .. .. .. ..@ plotOrder: int 1
## .. .. .. ..@ labpt : num [1:2] 356537 506466
## .. .. .. ..@ ID : chr "4"
## .. .. .. ..@ area : num 94710
## .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
## .. .. .. ..@ Polygons :List of 1
## .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
## .. .. .. .. .. .. ..@ labpt : num [1:2] 357005 506487
## .. .. .. .. .. .. ..@ area : num 184630
## .. .. .. .. .. .. ..@ hole : logi FALSE
## .. .. .. .. .. .. ..@ ringDir: int 1
## .. .. .. .. .. .. ..@ coords : num [1:7, 1:2] 356715 357153 357253 356768 356791 ...
## .. .. .. .. .. .. .. ..- attr(*, "dimnames")=List of 2
## .. .. .. .. .. .. .. .. ..$ : NULL
## .. .. .. .. .. .. .. .. ..$ : chr [1:2] "x" "y"
## .. .. .. ..@ plotOrder: int 1
## .. .. .. ..@ labpt : num [1:2] 357005 506487
## .. .. .. ..@ ID : chr "5"
## .. .. .. ..@ area : num 184630
## .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
## .. .. .. ..@ Polygons :List of 1
## .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
## .. .. .. .. .. .. ..@ labpt : num [1:2] 357513 506555
## .. .. .. .. .. .. ..@ area : num 99957
## .. .. .. .. .. .. ..@ hole : logi FALSE
## .. .. .. .. .. .. ..@ ringDir: int 1
## .. .. .. .. .. .. ..@ coords : num [1:6, 1:2] 357635 357640 357619 357406 357370 ...
## .. .. .. .. .. .. .. ..- attr(*, "dimnames")=List of 2
## .. .. .. .. .. .. .. .. ..$ : NULL
## .. .. .. .. .. .. .. .. ..$ : chr [1:2] "x" "y"
## .. .. .. ..@ plotOrder: int 1
## .. .. .. ..@ labpt : num [1:2] 357513 506555
## .. .. .. ..@ ID : chr "6"
## .. .. .. ..@ area : num 99957
## .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
## .. .. .. ..@ Polygons :List of 1
## .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
## .. .. .. .. .. .. ..@ labpt : num [1:2] 357146 505951
## .. .. .. .. .. .. ..@ area : num 326082
## .. .. .. .. .. .. ..@ hole : logi FALSE
## .. .. .. .. .. .. ..@ ringDir: int 1
## .. .. .. .. .. .. ..@ coords : num [1:6, 1:2] 356951 356573 357117 357586 357518 ...
## .. .. .. .. .. .. .. ..- attr(*, "dimnames")=List of 2
## .. .. .. .. .. .. .. .. ..$ : NULL
## .. .. .. .. .. .. .. .. ..$ : chr [1:2] "x" "y"
## .. .. .. ..@ plotOrder: int 1
## .. .. .. ..@ labpt : num [1:2] 357146 505951
## .. .. .. ..@ ID : chr "7"
## .. .. .. ..@ area : num 326082
## ..@ plotOrder : int [1:7] 7 5 2 3 6 1 4
## ..@ bbox : num [1:2, 1:2] 356279 505712 357640 507383
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : chr [1:2] "x" "y"
## .. .. ..$ : chr [1:2] "min" "max"
## ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
## .. .. ..@ projargs: chr "+proj=tmerc +lat_0=0 +lon_0=19 +k=0.9993 +x_0=500000 +y_0=-5300000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"plot(inter_p1_p2)
union_p1_p2 <- union(poligon1,poligon2) # unia z biblioteki raster
plot(union_p1_p2)
union_p1_p2 <- gUnion(poligon1,poligon2,unaryUnion_if_byid_false=FALSE) # złączeniem poligonów
plot(union_p1_p2)
diff_p1_p2 <- gDifference(poligon1,poligon2)
plot(diff_p1_p2)
diff_p2_p1 <- gDifference(poligon2,poligon1) #różnica nie jest operacją symetryczną
plot(diff_p2_p1)
xor_p1_p2 <- gSymdifference(poligon1,poligon2) #unia - przecięcie
plot(xor_p1_p2,col=2)
plot(poligon1)+plot(linie,add=TRUE) #usunięcie linii na przecięciu z poligonem
## numeric(0)xor_p1_linie <- gSymdifference(poligon1,linie)
plot(xor_p1_linie)
Raster
Raster jest regularną siatką zdefiniowanej zasięgiem, rozdzielczością i ilością wierszy i kolumn i w związku z tym mechanizmy selekcja może się opierać na tych samych zasadach jak w przypadku obiektów wektorowych tylko dla selekcji na podstawie wierszy i kolumn lub na podstawie zasięgu przestrzennego. Taka operacja nazywana jest przycinaniem, i pakiet raster posiada odpowiednie narzędzia, jak na przykład funkcję crop(). Dodatkowo, w przypadku zastosowania obiektów wielowarstwowych można zastosować procedurę selekcji poszczególnych warstw (bands) ze zbiorów wielowarstwowych, do czego służy funkcja subset(). W przypadku selekcji wybranych komórek istnieje kilka metod, pozwalających na operacje selekcji (wskazywania) które komórki spełniają wymagane kryterium.
Selekcja komórek na podstawie zasiegu przestrzennego
Jeżeli chcemy utworzyć nowy raster na podstawie zasięgu zdefiniowanego tylko przez numery wierszy i kolumn należy wykonać kilka kroków: zdefiniować extent; użyć funkcji crop(). Takie dziania są niezbędne, ze względu na konieczność uaktualnienia metadanych obiektu raster.
library(raster)
rast <- raster("DATA/pop_low.tif")
#selekcja na podstawie indeksów wierszy i kolumn
sel <- rast[10:20,50:70] #wybiera jedynie komórki, to nas teraz nie interesuje
ext <- extent(rast,105,121,78,99) #podajemy indeksy wieszy i kolumn
sel <- crop(rast,ext)
plot(sel)
ext## class : Extent
## xmin : 971565
## xmax : 976845
## ymin : 1840065
## ymax : 1844145#selekcja na podstawie zasięgu przestrzennego
ext <- extent(971565,976845,1840065,1844145)
sel <- crop(rast,ext)
plot(sel)
Selekcji komórek możemy też dokonać dowolnym innym obiektem geo-przestrzennym (klasy sp* i klasy raster). W przypadku obiektu klasy Sp do wykonania selekcji użyty zostanie ich bounding box.
#selekcja przez inny raster
selection <- rasterFromCells(rast,cellFromRowColCombine(rast,105:121,78:99)) # budujemy raster zgodny z wcześniej zdefiniowanym extent
sel <- intersect(rast,selection) #użyjemy ten raster bo nie mamy innego
plot(sel)
#selekcja przez obiekt wektorowy
selection <- readOGR("DATA/simple","lines")## OGR data source with driver: ESRI Shapefile
## Source: "DATA/simple", layer: "lines"
## with 1 features
## It has 1 fieldsplot(rast)+plot(selection,lwd=2,col=4,add=TRUE)
## numeric(0)sel <- intersect(rast,selection)
plot(sel)+plot(selection,lwd=2,col=4,add=TRUE)
## numeric(0)Selekcja na podstawie atrybutów (przez reklasyfikację do TRUE/FALSE)
W przypadku obiektu rastrowego pewnym rodzajem selekcji jest wskazanie, które komórki spełniają kryterium selekcji, a które nie i zapisanie wyniku jako warstwy logicznej. Do tego celu nalepiej wykorzystać funkcję calc(). W przypadku łączenia ze sobą operacji przestrzennych (a więc zmieniających geometrię) i atrybutowych (zwracających TRUE/FALSE)
sel <- rast[rast>200]
sel <- calc(rast,fun=function(x) {x>100}) # wykorzystujemy calc i funkcję operacji porównania, która zwraca prawdę i fałsz
summary(sel) # raster przechowuje wartości TRUE i FALSE ale wyświetla je jako 1 i 0## layer
## Min. 0
## 1st Qu. 0
## Median 0
## 3rd Qu. 0
## Max. 1
## NA's 6611v=getValues(sel) # jeżeli pobierzemy wartości, to zobaczymy że są typu TRUE/FALSE i mogą być wykorzystane jako wektor logiczny przy innnych operacjach
summary(v)## Mode FALSE TRUE NA's
## logical 16115 2801 6611Reklasyfikacja danych rastrowych
Re-klasyfikacja czyli zamiana oryginalnych wartości na klasy (kategorie) jest formą selekcji, która pokazuje zasięg przestrzenny obszarów spełniających pewne kryterium. Re-klasyfikację można wykonać na trzy sposoby: poprzez funkcję cut() tak jak dowolny zbiór wartości w R; funkcję reclasify() przy pomocy \(3 \times n\) macierzy re-klasyfikacji lub: w przypadku prostej zmiany wartości, przy pomocy funkcji subs().
Najprostsza forma re-klasyfikacji to użycie funkcji cut(), która podzieli oryginalny wektor danych zgodnie z liniami cięcia.
library(raster)
rast <- raster("DATA/pop_low.tif")
breaks <- c(minValue(rast),1,2,4,8,16,32,128,maxValue(rast))
rast_cut <- cut(rast,breaks=breaks,include.lowest=TRUE)
plot(rast_cut)
Funkcja subs pozwala użyć macierzy dwukolumnowej do zamiany jednych wartości na inne.
sub_mat <- data.frame(id=1:8,sub=c(1,2,3,4,5,6,10,20))
rast_sub <- subs(rast_cut,sub_mat)
plot(rast_sub)
Macierz re-klasyfikacji nie musi być w postaci macierzy: wystarczy że powtarza się sekwencja: pierwsza wartość oznacza dolny zakres przedziału; druga – górny zakres przedziału a trzecia wartość wartość do której następuje re-klasyfikacja. Re-klasyfikacja tworzy nową warstwę numeryczną (a nie integer). Aby re-klasyfikacja była udana, należy rozpocząć ją od najmniejszej wartości a skończyć na największej: przy pomocy funkcji minValue()-1 i maxValue()+1.
rec_mat <- c(minValue(rast)-1,1,1,
1,2,2,
2,4,3,
4,8,4,
8,16,5,
16,32,6,
32,128,10,
128,maxValue(rast)+1,20)
rast_rec <- reclassify(rast,rec_mat)
plot(rast_rec)