Wstęp

Podstawowe obiekty geo-przestrzenne to punkty, linie i poligony oraz siatki rastrowe (grid). Do obsługi obiektów geo-przestrzennych istnieje kilka bibliotek, z których trzy najważniejsze to:

Wymagane biblioteki: rgdal, sp, raster Wymagana wiedza: podstawy programowania, systemy informacji przestrzennej

Tworzenie obiektów geoprzestrzennych z istniejących zbiorów danych geoprzestrzennych

Czytanie do obiektów klasy Spatial*

To czytania/tworzenia danych geo-przestrzennych klasy Spatial* służy biblioteka rgdal. Obiekty tworzymy przy pomocy funkcji readOGR (dane wktorowe) lub readGDAL (dane rastrowe). Biblioteka rgdal pozwala wczytać do pamięci dowolny plik geoprzestrzenny obsługiwany przez bibliotekę GDAL/OGR. Wynikiem jest obiekt klasy _Spatial*DataFrame_, przechowywany w pamięci operacyjnej. Z tego powodu należy w pierwszej kolejności sprawdzić, czy wczytywany zbiór danych nie jest zbyt duży. Funkcja GDALInfo zwraca wiele informacji na temat zbioru danych w tym również liczbę wierszy i kolumn.

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
punkty <- readOGR(dsn="DATA/simple", layer="points") #czytanie obiektów wektorowych
## OGR data source with driver: ESRI Shapefile 
## Source: "DATA/simple", layer: "points"
## with 9 features
## It has 2 fields
linie <- readOGR(dsn="DATA/simple", layer="lines")
## OGR data source with driver: ESRI Shapefile 
## Source: "DATA/simple", layer: "lines"
## with 1 features
## It has 1 fields
meta <- GDALinfo("DATA/pop_low.tif") #metadane obiektu rastrowego
meta["rows"]*meta["columns"] #rozmiar pliku wiersze * columna
##  rows 
## 25527
rast <- readGDAL("DATA/pop_low.tif")
## DATA/pop_low.tif has GDAL driver GTiff 
## and has 127 rows and 201 columns
image(rast,col=rainbow(20)) #wyświetlenie obiektu rastrowego
points(punkty) #wyświetlenie wektorów
lines(linie)

Struktura wewnętrzna obiektów geoprzestrzennych klasy Spatial*

Obiekty wektorowe klasy spatial*DataFrame składają się z czterech głównych slotów:

  • @data – przechowuje data.frame o liczbie wierszy równej ilości obiektów geometrycznych
  • @geometry – przechowuje tablicę punktów lub listę linii albo poligonów.
  • @bbox – przechowuje zasięg (extent) danych geo-przestrzennych. Najczęściej wyznaczany maksymalne/minimalne wartości współrzędnych
  • @crs – przechowuje informację o układzie geo-przestrzennym
str(punkty)
## Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
##   ..@ data       :'data.frame':  9 obs. of  2 variables:
##   .. ..$ id : int [1:9] 1 2 3 4 5 6 7 8 9
##   .. ..$ id2: int [1:9] 2 4 6 8 10 12 14 16 18
##   ..@ coords.nrs : num(0) 
##   ..@ coords     : num [1:9, 1:2] 956936 958663 968119 979672 976753 ...
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : NULL
##   .. .. ..$ : chr [1:2] "coords.x1" "coords.x2"
##   ..@ bbox       : num [1:2, 1:2] 956936 1846136 996117 1863321
##   .. ..- 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=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"| __truncated__

Struktura obiektu SpatialGridDataFrame jest bardzo podobna do struktury innych obiektów typu Spatial*. Posiada 4 sloty:

  • @data – przechowuje data.frame o liczbie wierszy równej ilości wczytanych pasm (bands)
  • @grid – przechowuje topologię siatki lub w przypadku SpatialPixelsDataframe #uzupełnić
  • @bbox – przechowuje zasięg (extent) danych geo-przestrzennych. Najczęściej wyznaczany maksymalne/minimalne wartości współrzędnych lub siatki
  • @crs – przechowuje informację o układzie geo-przestrzennym
str(rast)
## Formal class 'SpatialGridDataFrame' [package "sp"] with 4 slots
##   ..@ data       :'data.frame':  25527 obs. of  1 variable:
##   .. ..$ band1: num [1:25527] NA NA NA NA NA NA NA NA NA NA ...
##   ..@ grid       :Formal class 'GridTopology' [package "sp"] with 3 slots
##   .. .. ..@ cellcentre.offset: Named num [1:2] 953205 1838745
##   .. .. .. ..- attr(*, "names")= chr [1:2] "x" "y"
##   .. .. ..@ cellsize         : num [1:2] 240 240
##   .. .. ..@ cells.dim        : int [1:2] 201 127
##   ..@ bbox       : num [1:2, 1:2] 953085 1838625 1001325 1869105
##   .. ..- 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=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"

Czytanie do obiektów klasy raster*

Unikalną zaletą obiektów klasy raster* jest możliwość przechowywania danych na dysku zamiast w pamięci operacyjnej. Ponieważ dane rastrowe często swoim rozmiarem przekraczają dostępną pamięć operacyjną brak możliwości obsługi dużych zbirów danych jest jednym z ograniczeń biblioteki sp. Biblioteka raster obsługuje następujące typy obiektów rastrowych:

  • RasterLayer – obiekt typu RasterLayer zawiera tylko jedną warstwę rastrową. Utworzymy przy pomocy funkcji raster()
  • RasterStack – służy do przechowywania danych rastrowych pochodzących z różnych plików ale mających ten sam początek i rozdzielczość. Obiekty RasteStack tworzymy poleceniem stack()
  • RasterBrick – pozwala na wczytywanie wielowarstwowych (multiband) plików rastrowych. Działa znacząco szybciej niż RasterStack. Obiekt typu RasterBrick tworzymy poleceniem brick.

W zależności od typu przetwarzanych danych należy wybrać typ obiektu docelowego. W przypadku danych jednowarstwowych typ RasterLayer i uniwersalna funkcja raster() jest najbardziej oczywistym wyborem. W przypadku, gdy chcemy przetwarzać dane wielowarstwowe wybór zależy od rodzaju danych wejściowych. Jeżeli dane znajdują się w osobnych plikach rozwiązaniem jest stack() natomiast dla plików wielowarstwowych (na przykład zdjęć multispektralych) polecenie odpowiednie jest polecenie brick(). Późniejsza konwersja pomiędzy typami obiektów jest jak najbardziej możliwa.

require(raster)
## Loading required package: raster
layer <- raster("DATA/pop_low.tif")
layer
## class       : RasterLayer 
## dimensions  : 127, 201, 25527  (nrow, ncol, ncell)
## resolution  : 240, 240  (x, y)
## extent      : 953085, 1001325, 1838625, 1869105  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
## data source : /home/jarekj/Dropbox/dydaktyka/prog_geo/PUBLISH/DATA/pop_low.tif 
## names       : pop_low 
## values      : 0, 1507.513  (min, max)

Struktura obiektu typu Raster* jest bardzo złożona, ze względu na konieczność zapewnienia obsługi danych na dysku, w tym informacji o strukturze zewnętrznego zbioru danych takiej jak typ danych, kolejność bitów, układ danych w pliku itp. Obiekty tworzone na podstawie istniejących danych geo-przestrzennych są domyślnie wczytywane do pamięci.

str(layer)
## Formal class 'RasterLayer' [package "raster"] with 12 slots
##   ..@ file    :Formal class '.RasterFile' [package "raster"] with 13 slots
##   .. .. ..@ name        : chr "/home/jarekj/Dropbox/dydaktyka/prog_geo/PUBLISH/DATA/pop_low.tif"
##   .. .. ..@ datanotation: chr "FLT8S"
##   .. .. ..@ byteorder   : chr "little"
##   .. .. ..@ nodatavalue : num -Inf
##   .. .. ..@ NAchanged   : logi FALSE
##   .. .. ..@ nbands      : int 1
##   .. .. ..@ bandorder   : chr "BIL"
##   .. .. ..@ offset      : int 0
##   .. .. ..@ toptobottom : logi TRUE
##   .. .. ..@ blockrows   : int 5
##   .. .. ..@ blockcols   : int 201
##   .. .. ..@ driver      : chr "gdal"
##   .. .. ..@ open        : logi FALSE
##   ..@ data    :Formal class '.SingleLayerData' [package "raster"] with 13 slots
##   .. .. ..@ values    : logi(0) 
##   .. .. ..@ offset    : num 0
##   .. .. ..@ gain      : num 1
##   .. .. ..@ inmemory  : logi FALSE
##   .. .. ..@ fromdisk  : logi TRUE
##   .. .. ..@ isfactor  : logi FALSE
##   .. .. ..@ attributes: list()
##   .. .. ..@ haveminmax: logi TRUE
##   .. .. ..@ min       : num 0
##   .. .. ..@ max       : num 1508
##   .. .. ..@ band      : int 1
##   .. .. ..@ unit      : chr ""
##   .. .. ..@ names     : chr "pop_low"
##   ..@ legend  :Formal class '.RasterLegend' [package "raster"] with 5 slots
##   .. .. ..@ type      : chr(0) 
##   .. .. ..@ values    : logi(0) 
##   .. .. ..@ color     : logi(0) 
##   .. .. ..@ names     : logi(0) 
##   .. .. ..@ colortable: logi(0) 
##   ..@ title   : chr(0) 
##   ..@ extent  :Formal class 'Extent' [package "raster"] with 4 slots
##   .. .. ..@ xmin: num 953085
##   .. .. ..@ xmax: num 1e+06
##   .. .. ..@ ymin: num 1838625
##   .. .. ..@ ymax: num 1869105
##   ..@ rotated : logi FALSE
##   ..@ rotation:Formal class '.Rotation' [package "raster"] with 2 slots
##   .. .. ..@ geotrans: num(0) 
##   .. .. ..@ transfun:function ()  
##   ..@ ncols   : int 201
##   ..@ nrows   : int 127
##   ..@ crs     :Formal class 'CRS' [package "sp"] with 1 slot
##   .. .. ..@ projargs: chr "+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
##   ..@ history : list()
##   ..@ z       : list()
inMemory(layer)
## [1] FALSE
layer <- readAll(layer)
inMemory(layer)
## [1] TRUE
layer
## class       : RasterLayer 
## dimensions  : 127, 201, 25527  (nrow, ncol, ncell)
## resolution  : 240, 240  (x, y)
## extent      : 953085, 1001325, 1838625, 1869105  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
## data source : in memory
## names       : pop_low 
## values      : 0, 1507.513  (min, max)

Jeżeli z jakiś względów zależy nam aby dane w całości były przechowywane w pamięci należy użyć funkcji readAll(). Natomiast dowolna konwersja między typami raster*, powoduje, że wynik jest zapisywany w pamięci.

Tworzenie obiektów geoprzestrzennych od podstaw

SpatialPoints i PointsDataFrame

Przy użyciu narzędzi z biblioteki sp utworzenie obiektu punktowego z danych geometrycznych jest bardzo proste, ponieważ istnieje tylko jeden obiekt: SpatialPoints, składający się z punktów. SpatialPointsDataframe to punkty z układem odniesienia. Potrzebne są: geometria w postaci macierzy dwóch kolumn, każda zawierająca współrzędne jednego punktu (na przykład losowo wygenerowane punkty), data.frame z atrybutami (na przykład sekwencja) oraz system referencji przestrzennej (nie jest wymagany ale będą problemy z projekcją danych). Bounding Box jeżeli nie został podany jest generowany automatycznie.

Tworzenie obiektów geo-przestrzennych od podstaw jest żmudną czynnością, polegającą na ręcznym wypisywaniu współrzędnych. Dlatego te obiekty są tworzone na podstawie obliczeń

library(sp)
punkty <- cbind(runif(20,520000,521000),runif(20,440000,444000))
dane <- data.frame(id=seq(1:20))
ref <- CRS("+init=epsg:2180")
sp_punkty <- SpatialPointsDataFrame(punkty,dane,proj4string=ref)
par(mfrow=c(1,2))
plot(punkty,asp=1)
plot(sp_punkty)

par(mfrow=c(1,1))

Jeżeli chcemy utworzyć tylko geometrię, (bez atrybutów) należy użyć polecenia SpatialPoints(). Danych tego typu nie będzie można wyeksportować przy użyciu funkcji writeOGR()

SpatialLines i SpatialLinesDataFrame

Przy użyciu biblioteki sp tworzenie linii jest dużo bardziej skomplikowane, gdyż występuje hierarchia obiektów (uwaga na wielkość liter:

  • Line – pojedyncza linia
  • Lines – grupa (lista) linii.
  • SpatialLines – Linie w układzie geo-przestrzennych

Podobnie jak w przypadku puntów SpatalLinesDataFrame to linie z układem odniesienia. Aby utworzyć SpatialLineDataframe, w pierwszej kolejności należy utworzyć geometrię. utworzymy dwie linie:

simple_line1 <- cbind(runif(20,520000,520200),sort(runif(20,440000,444200))) #sort ma na celu wygenerowanie kolejności
simple_line2 <- cbind(sort(runif(20,519000,521000)),sort(runif(20,440000,444100)))
plot(simple_line1,type="l",asp=1)
lines(simple_line2,col=2)

Następny krok to promocja listy punktów do obiektu typu Line, a następnie obiektu Lines. Zależy nam aby każda linia miała osobny identyfikator.

line_geom1 <- Line(simple_line1)
line_geom2 <- Line(simple_line2)
str(line_geom1)
## Formal class 'Line' [package "sp"] with 1 slot
##   ..@ coords: num [1:20, 1:2] 520156 520106 520159 520088 520118 ...
#Promocja do klasy Lines i nadanie Identyfikatora obiektu
lines_geom1 <- Lines(list(line_geom1),ID="1")
lines_geom2 <- Lines(list(line_geom2),ID="2")
str(lines_geom1)
## Formal class 'Lines' [package "sp"] with 2 slots
##   ..@ Lines:List of 1
##   .. ..$ :Formal class 'Line' [package "sp"] with 1 slot
##   .. .. .. ..@ coords: num [1:20, 1:2] 520156 520106 520159 520088 520118 ...
##   ..@ ID   : chr "1"

To w sumie nieco dziwne postępowanie wynika, z faktu że pojedyncza linia może nie być wcale pojedynczym obiektem geo-przestrzennym (a więc posiadającym indywidualne atrybuty). Poprzerywane poziomice są dobrym przykładem obiektu przestrzennego składającego się z kilku geometrii ale mającym wspólny atrybut (rzędną). Konieczność budowania zagnieżdżonych list jest efektem natury języka R. Złożone obiekty mogą być tworzone tylko jako listy.

Mając obiekty klasy Lines, możemy zamienić je na obiekt klasy SpatialLines – poprzez nadanie uprzednio utworzonego CRS a następnie do klasy SpatialLinesDataFrame, poprzez przypisanie tabeli atrybutów. Tabela atrybutów powinna liczyć tyle elementów, ile mamy osobnych obiektów geo-przestrzennych (w własnym ID) a więc dwóch.

SpatialLines_geom <- SpatialLines(list(lines_geom1,lines_geom2),proj4string=ref)
data <- data.frame(id=c(1,2),color=c("red","black"))
sp_lines <- SpatialLinesDataFrame(SpatialLines_geom,data)
plot(sp_lines)

str(sp_lines)
## Formal class 'SpatialLinesDataFrame' [package "sp"] with 4 slots
##   ..@ data       :'data.frame':  2 obs. of  2 variables:
##   .. ..$ id   : num [1:2] 1 2
##   .. ..$ color: Factor w/ 2 levels "black","red": 2 1
##   ..@ lines      :List of 2
##   .. ..$ :Formal class 'Lines' [package "sp"] with 2 slots
##   .. .. .. ..@ Lines:List of 1
##   .. .. .. .. ..$ :Formal class 'Line' [package "sp"] with 1 slot
##   .. .. .. .. .. .. ..@ coords: num [1:20, 1:2] 520156 520106 520159 520088 520118 ...
##   .. .. .. ..@ ID   : chr "1"
##   .. ..$ :Formal class 'Lines' [package "sp"] with 2 slots
##   .. .. .. ..@ Lines:List of 1
##   .. .. .. .. ..$ :Formal class 'Line' [package "sp"] with 1 slot
##   .. .. .. .. .. .. ..@ coords: num [1:20, 1:2] 519115 519219 519287 519319 519385 ...
##   .. .. .. ..@ ID   : chr "2"
##   ..@ bbox       : num [1:2, 1:2] 519115 440262 520966 444192
##   .. ..- 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 "+init=epsg:2180 +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"| __truncated__

Pakiet raster zawiera uproszczony mechanizm budowy obiektów linijnych, przy pomocy funkcji spLines(). Uproszczenie zakłada że każda geometria to osobny obiekt. Identyfikator jest przypisywany automatycznie. Dane muszą być zgodne z liczbą obiektów:

library(raster)
sp_lines_rast <- spLines(simple_line1,simple_line2,attr=data,crs=ref)
plot(sp_lines_rast,col=c(1,2))

W przypadku, gdy pojedyncze obiekty geometryczne chcemy grupować, należy w sposób jawny utworzyć z nich listy:

simple_line3 <- cbind(runif(20,522000,522200),sort(runif(20,440000,444200))) #dodatkowa linia
sp_lines_rast2 <- spLines(list(simple_line1,simple_line2),simple_line3,attr=data,crs=ref) #ale nadal mamy tylko dwa obiekty więc data frame ma dwa atrybuty
plot(sp_lines_rast2,col=c(1,2))

SpatialPolygons i SpatialPolygonsDataFrame

Struktura wewnętrzna poligonów jest bardzo zbliżona, do linii. Jedyną różnicą jest to ze mogą one zawierać dziury. Obiekty wewnętrzne poligonów to:

  • Polygon – pojedynczy pierścień (zamknięty ciąg puntów)
  • Polygons – Grupa pierścieni, z reguły pierścień główny i dziury.
  • SpatialPolgons – Poligony w układzie geo-przestrzennym

Podobnie jak w przypadku linii SpatalPolygonsDataFrame to Poligony z układem odniesienia. Aby utworzyć Poligony ożyjemy metody uproszczonej (raster::spPolygons). Utworzymy trzy poligony, pojedynczy, z otworem i z wyspą.

simple_polygon1 <- rbind(c(0,50),c(10,60),c(30,70),c(70,50),c(40,55))

simple_polygon2 <- rbind(c(40,0),c(40,40),c(80,40),c(80,0))
hole2 <- rbind(c(50,10),c(70,10),c(70,30),c(50,30)) #dziury muszą mieć kolejność punktów przeciwnie do ruchy wskazówek zegara

simple_polygon3 <- rbind(c(0,20),c(20,30),c(35,20),c(20,10))
island3 <- rbind(c(0,5),c(10,10),c(5,5))

data <- data.frame(ID=c(1,2,3),type=c("ring","dziura","wyspa"))

sp_polygon_rast <- spPolygons(simple_polygon1, list(simple_polygon2,hole2), 
                              list(simple_polygon3,island3),
                              attr=data,crs=ref)
plot(sp_polygon_rast,col=c(1,2,4))

Alternatywnie do definiowania dziur w poligonach można wymusić typ poligonu używając funkcji Polygon bezpośrednio. Zagwarantuje to nam odpowiednią kolejność punktów.

hole2 <- Polygon(hole2,hole=TRUE)

SpatialGrid i SpatialGridDataFrame

Rastrowy obiekt geo-przestrzenny klasy SpatialGrid można utworzyć na dwa sposoby: poprzez zdefiniowanie topologii siatki (a więc geometrii) lub poprzez przekształcenia zbioru punktów na strukturę SpatialPixels. To ostatnie rozwiązanie stosujemy do niekompletnych siatek, zdefiniowanych siatką regularnych punktów.

Tworzenie topologii.

W celu zdefiniowanie topologii musimy znać trzy parametry: Początek siatki, wielkość komórki, rozmiar siatki. Należy przy tym pamiętać że początek siatki jest przesunięty o \(1\over2\) wielkości komórki. Dla załączonego powyżej przykładu będą to następujące parametry:

topology <- GridTopology(
  cellcentre.offset = c(953205,1838745),
  cellsize=c(240,240),
  cells.dim=c(201,127)
)
siatka <- SpatialGrid(topology,proj4string = CRS("+init=epsg:5070"))

W następnym kroku możemy dodać do istniejącej topologii dowolne dane. Dane muszą być w formacie data.frame a ich wielkość musi być zgodna z rozmiarem siatki. Dodamy trzy warstwy (bands): pustą, kolejkę oraz dane losowe od 0 do 1. W ten sposób zostanie utworzony gotowy obiekt geo-przestrzenny.

require(sp)
size <- siatka@grid@cells.dim[1]*siatka@grid@cells.dim[2]
size
## [1] 25527
dane <- data.frame(empty=rep(NA,size),sequ=seq(1:size),los=runif(size))
siatka <- SpatialGridDataFrame(siatka,dane)
par(mfcol=c(1,2))
image(siatka,attr=2)
image(siatka,attr=3) 

par(mfcol=c(1,1))
str(siatka)
## Formal class 'SpatialGridDataFrame' [package "sp"] with 4 slots
##   ..@ data       :'data.frame':  25527 obs. of  3 variables:
##   .. ..$ empty: logi [1:25527] NA NA NA NA NA NA ...
##   .. ..$ sequ : int [1:25527] 1 2 3 4 5 6 7 8 9 10 ...
##   .. ..$ los  : num [1:25527] 0.934 0.118 0.126 0.384 0.582 ...
##   ..@ grid       :Formal class 'GridTopology' [package "sp"] with 3 slots
##   .. .. ..@ cellcentre.offset: num [1:2] 953205 1838745
##   .. .. ..@ cellsize         : num [1:2] 240 240
##   .. .. ..@ cells.dim        : int [1:2] 201 127
##   ..@ bbox       : num [1:2, 1:2] 953085 1838625 1001325 1869105
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : NULL
##   .. .. ..$ : chr [1:2] "min" "max"
##   ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
##   .. .. ..@ projargs: chr "+init=epsg:5070 +proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS8"| __truncated__

SpatialPixels. Tworzenie siatki na podstawie punktów

Jeżeli mamy zestaw regularnych punktów możemy utworzyć obiekt klasy SpatialPixels lub SpatialPixelsDataFrame. Do zbudowania obiektu SpatialPixel* wykozystamy macierz. Aby utworzyć SpatialPixelDataFrame w piwrwszej kolejności należy utworzyć obiekt klasy SpatialPoint*

m <- matrix(seq(1:400),20,20)
m[1:5,1:3] <- NA #wycinamy prostokąt
m[17:20,14:20] <- NA
m[m<300] <- 1
m[m>=300] <- 2
image(m)

g <- expand.grid(m,xc=1:20,yc=1:20)
str(g)
## 'data.frame':    160000 obs. of  3 variables:
##  $ Var1: num  NA NA NA NA NA 1 1 1 1 1 ...
##  $ xc  : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ yc  : int  1 1 1 1 1 1 1 1 1 1 ...
##  - attr(*, "out.attrs")=List of 2
##   ..$ dim     : Named int  400 20 20
##   .. ..- attr(*, "names")= chr  "" "xc" "yc"
##   ..$ dimnames:List of 3
##   .. ..$ Var1: chr  "Var1=NA" "Var1=NA" "Var1=NA" "Var1=NA" ...
##   .. ..$ xc  : chr  "xc= 1" "xc= 2" "xc= 3" "xc= 4" ...
##   .. ..$ yc  : chr  "yc= 1" "yc= 2" "yc= 3" "yc= 4" ...
coordinates(g) = ~xc+yc #ależy stosować znak równości
summary(g)
## Object of class SpatialPointsDataFrame
## Coordinates:
##    min max
## xc   1  20
## yc   1  20
## Is projected: NA 
## proj4string : [NA]
## Number of points: 160000
## Data attributes:
##       Var1      
##  Min.   :1.000  
##  1st Qu.:1.000  
##  Median :1.000  
##  Mean   :1.224  
##  3rd Qu.:1.000  
##  Max.   :2.000  
##  NA's   :17200
gridded(g) <- TRUE
summary(g)
## Object of class SpatialPixelsDataFrame
## Coordinates:
##    min  max
## xc 0.5 20.5
## yc 0.5 20.5
## Is projected: NA 
## proj4string : [NA]
## Number of points: 160000
## Grid attributes:
##    cellcentre.offset cellsize cells.dim
## xc                 1        1        20
## yc                 1        1        20
## Data attributes:
##       Var1      
##  Min.   :1.000  
##  1st Qu.:1.000  
##  Median :1.000  
##  Mean   :1.224  
##  3rd Qu.:1.000  
##  Max.   :2.000  
##  NA's   :17200

Tworzenie obiektów klasy Raster* od podstaw

Funkcja raster() (oraz stack() i brick) jest funkcją polimorficzną tzn. może wykonywać wiele czynności, w zależności od typów i rodzaju przekazanych danych. W przypadku jeżeli jest to ścieżka dostępu do istniejącego zbioru danych geo-przestrzennych funkcja po prostu utworzy obiekt z danymi przechowywanymi na dysku. Jeżeli chcemy utworzyć obiekt klasy Raster* od podstaw musimy podać podstawowe informacje na temat planowanego obiektu. Można też utworzyć domyślny obiekt który będzie obejmował całą kulę ziemską, w układzie WGS84 o rozdzielczości 1 stopnia:

raster()
## class       : RasterLayer 
## dimensions  : 180, 360, 64800  (nrow, ncol, ncell)
## resolution  : 1, 1  (x, y)
## extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0

Tworząc dane o innym zasięgu przestrzennym, musimy znać zakres siatki (bez przesunięcia jak w przypadku sp) oraz rozdzielczość lub ilość kolumn i wierszy. Brak georferencji spowoduje przyjęcie domyślnej WGS84. Z tego powodu wartości topologii siatki muszą być wyliczone przed przesłaniem ich do programu. Jeżeli nie chcemy podawać, georeferencji parametr crs należy ustawić na NA.

xmin <- 500000
ymin <- 400000
nrow <- 100
ncol <- 100
res <- 10
ext <- extent(c(xmin,xmin+ncol*res,ymin,ymin+nrow*res)) #utworzenie obiektu extent
rast <- raster(ext,ncol=ncol,nrow=nrow,res=res,crs=NA)
rast
## class       : RasterLayer 
## dimensions  : 100, 100, 10000  (nrow, ncol, ncell)
## resolution  : 10, 10  (x, y)
## extent      : 5e+05, 501000, 4e+05, 401000  (xmin, xmax, ymin, ymax)
## coord. ref. : NA

Dowolne parametry geometryczne obiektu typu Raster* mogą być zmienione poprzez zmianę wartości atrybutów. Funkcja automatycznie przeliczy pozostałe elementy topologii. Wszelkie zmiany w topologii należy wykonać przed dodaniem wartości. Późniejsza ich zmiana może spowodować usunięcie warstwy danych.

rast
## class       : RasterLayer 
## dimensions  : 100, 100, 10000  (nrow, ncol, ncell)
## resolution  : 10, 10  (x, y)
## extent      : 5e+05, 501000, 4e+05, 401000  (xmin, xmax, ymin, ymax)
## coord. ref. : NA
nrow(rast)
## [1] 100
nrow(rast) <- 50
rast
## class       : RasterLayer 
## dimensions  : 50, 100, 5000  (nrow, ncol, ncell)
## resolution  : 10, 20  (x, y)
## extent      : 5e+05, 501000, 4e+05, 401000  (xmin, xmax, ymin, ymax)
## coord. ref. : NA
res(rast) <- 10
rast
## class       : RasterLayer 
## dimensions  : 100, 100, 10000  (nrow, ncol, ncell)
## resolution  : 10, 10  (x, y)
## extent      : 5e+05, 501000, 4e+05, 401000  (xmin, xmax, ymin, ymax)
## coord. ref. : NA

Aby dodać wartości należy wykonać funkcję values(). Funkcja ta zwraca wartości lub powoduje przypisanie wartości do obiektu. Wartości można zobaczyć wykonując funkcję plot() lub image. Należy przypomnieć, że w pakiecie sp do rysowania mapy służy funkcja funkcjaimage() lub w ograniczonym zakresie spplot().

head(values(rast)) # funkcja head zkraca wiekość danych wyjściowych
## [1] NA NA NA NA NA NA
values(rast) <- runif(ncell(rast))
image(rast)

Konwerjsja z- i do- innych typów obietów geoprzestrzennych do obiektu Raster*

Obiekty klasy Raster* mogą być tworzone na podstawie innych obiektów geoprzestrzennych oraz macierzy: SpatialGrid, SpatialPixel im (pakiet spatstat) a także z macierzy. W przypadku tworzenia obiektu raster z macierzy należy podać zasięg geo-przestrzenny obiektu. W tym przykładzie wykorzystamy utworzony wcześniej obiekt ext. W związku z tym nasza macierz powinna mieć rozmiar \(100\times100\). Ilość kolumn i wierszy zostanie pobrana z macierzy, natomiast rozdzielczość zostanie wyliczona automatycznie.

m <- matrix(rnorm(10000),100,100)
rast <- raster(m)
extent(rast) <- ext
rast
## class       : RasterLayer 
## dimensions  : 100, 100, 10000  (nrow, ncol, ncell)
## resolution  : 10, 10  (x, y)
## extent      : 5e+05, 501000, 4e+05, 401000  (xmin, xmax, ymin, ymax)
## coord. ref. : NA 
## data source : in memory
## names       : layer 
## values      : -4.230765, 3.602173  (min, max)
par(mfcol=c(1,2))
image(m)
image(rast)

par(mfcol=c(1,1))

Do zilustrowania konwersji z obiektu SpatialGrid* użyjemy wcześniej utworzonego obiektu siatka. W przypadku gdy obiekt SpatialGridDataFrame posiada więcej niż jedną warstwę danych należy wskazać, która warstwa ma być użyta w obiekcie RasterLayer lub wybrać RasterBrick. W naszym przykładzie w obiekcie siatka pierwsza warstwa jest pusta, pozostałe to warstwa losowa i sekwencja. Dodatkowo jeżeli chcemy aby wynik konwersji był przechowywany na dysku (ze względu na rozmiar) należy zastosować parametr filename. Parametr filename nie może być stosowany w każdej konwersji. Szczegółowo kwestie związane z zarządzaniem pamięcią zostaną omówione w rozdziale poświęconym tworzeniu funkcji przetwarzających dane rastrowe.

rast_grid <- raster(siatka)
plot(rast_grid)

rast_grid <- raster(siatka,layer=3) #warstwa losowa
plot(rast_grid)

brick_grid <- brick(siatka[2:3]) # pomijamy pierwszą warstwę pustą
plot(brick_grid)

brick_grid # in memory
## class       : RasterBrick 
## dimensions  : 127, 201, 25527, 2  (nrow, ncol, ncell, nlayers)
## resolution  : 240, 240  (x, y)
## extent      : 953085, 1001325, 1838625, 1869105  (xmin, xmax, ymin, ymax)
## coord. ref. : +init=epsg:5070 +proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0 
## data source : in memory
## names       :         sequ,          los 
## min values  : 1.000000e+00, 6.180955e-06 
## max values  : 2.552700e+04, 9.999231e-01
brick_grid <- brick(raster(siatka,layer=2),raster(siatka,layer=3),filename="DATA/tmp3.tif",overwrite=TRUE)
brick_grid # on disk
## class       : RasterBrick 
## dimensions  : 127, 201, 25527, 2  (nrow, ncol, ncell, nlayers)
## resolution  : 240, 240  (x, y)
## extent      : 953085, 1001325, 1838625, 1869105  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
## data source : /home/jarekj/Dropbox/dydaktyka/prog_geo/PUBLISH/DATA/tmp3.tif 
## names       :       tmp3.1,       tmp3.2 
## min values  : 1.000000e+00, 6.180955e-06 
## max values  : 2.552700e+04, 9.999231e-01

Biblioteka raster posiada narzędzia do bezpośredniej konwersji obiektów z biblioteki spatstat, służącej do analiz gęstościowych punktów oraz kilku innych mniej popularnych specjalistycznych narzędzi. Do konwersji na obiekty klasy Spatial* służy funkcja as(). Podobnie funkcja as() z biblioteki maptools. pozwala na konwersję pomiędzy obiektami spatal* i spatstat Ten mechanizm zostanie omówiony przy okazji omawiania analiz gęstościowych.

Zapis obiektów geoprzestrzennych

Do zapisu obiektów wektorowych klasy spatial* używamy funkcji writeOGR(). Funkcja ta wymaga zdefiniowania, żródła danych. Na tym etapie będzie to katalog, gdzie zapiszemy obkiet typu ESRI Shapefile. Przy pomocy funkcji writeOGR() możemy zapisać tylko obiekty SpatialDataFrame, nie możemy zapisać obiektów pozbawianych warstwy danych (a więc zwykłych Spatial)

dsn <- "DATA" #katalog data, ścieżka względna
driver="ESRI Shapefile"
writeOGR(sp_lines,dsn=dsn,layer="sp_linie",driver=driver,overwrite_layer = TRUE) #nie dajemy rozszerzenia
writeOGR(sp_polygon_rast,dsn=dsn,layer="sp_polygon",driver=driver,overwrite_layer = TRUE)

Do zapisu obiektów klasy Raster* służy funkcja writeRaster() z biblioteki raster natomiast obiekty SpatialGridDataFrame zapisywane są funkcją writeGDAL() z biblioteki rgdal. Obiektów Raster* przechowywanych na dysku nie trzeba zapisywać. Należy pamiętać, że funkcja writeGDAL() może nadpisać dane. Funkcja writeRaster wymaga ustawienia flagi overwrite na TRUE

writeRaster(rast_grid,filename="DATA/rast_grid.tif",overwrite=TRUE)
writeGDAL(siatka[2:3],fname="DATA/sp_siatka.tif") #nadpisuje istniejące dane