2 R a dane przestrzenne

2.1 Wprowadzenie

2.1.1 Reprezentacja danych nieprzestrzennych

  • Wektory (ang. vector):
    • liczbowe (ang. integer, numeric) - c(1, 2, 3) i c(1.21, 3.32, 4.43)
    • znakowe (ang. character) - c("jeden", "dwa", "trzy")
    • logiczne (ang. logical) - c(TRUE, FALSE)
    • czynnikowe (ang. factor) - c("jeden", "dwa", "trzy", "jeden")
  • Ramki danych (ang. data.frame) - to zbiór zmiennych (kolumn) oraz obserwacji (wierszy) zawierających różne typy danych
  • Macierze (ang. matrix)
  • Listy (ang. list)

2.1.2 Pakiety

R zawiera wiele funkcji pozwalających na przetwarzanie, wizualizację i analizowanie danych przestrzennych. Zawarte są one w szeregu pakietów (zbiorów funkcji), między innymi:

  • GIS - sf, raster, sp, rgdal, rgeos
  • Geostatystyka - gstat, geoR, geoRglm

Więcej szczegółów na ten temat pakietów R służących do analizy przestrzennej można znaleźć pod adresem https://cran.r-project.org/web/views/Spatial.html.

2.1.3 Reprezentacja danych przestrzennych

Dane przestrzenne mogą być reprezentowane w R poprzez szereg różnych klas obiektów z użyciem różnych pakietów R. Przykładowo dane wektorowe mogą być reprezentowane poprzez obiekty klas Spatial* z pakietu sp oraz obiekty klasy sf z pakietu sf.

Wszystkie obiekty klasy Spatial* z pakietu sp zawierają tablę atrybutów oraz dwie dodatkowe informacje:

- bounding box (`bbox`) - obwiednia - określa zasięg danych
- CRS (`proj4string`) - układ współrzędnych

Ten pakiet definiuje klasę obiektów - sposób reprezentacji danych. Najczęściej stosowane obiekty klasy Spatial* to SpatialPointsDataFrame, SpatialPolygonsDataFrame oraz SpatialGridDataFrame. Ostatnia klasa reprezentuje dane rastrowe.
Dodatkowo ten pakiet współpracuje z pakietami rgdal służącym do wczytywania i zapisywania danych oraz rgeos służącym do przetwarzania danych przestrzennych. W oparciu o pakiet sp powstało kilkaset dodatkowych pakietów R do analizy danych przestrzennych.

Dane rastrowe są reprezentowane między innymi poprzez klasę SpatialGridDataFrame z pakietu sp oraz obiekty klasy Raster* z pakietu raster, tj. RasterLayer, RasterStack, RasterBrick.

Więcej o pakiecie sf oraz raster można przeczytać w rozdziale drugim Geocomputation with R.

2.1.4 GDAL/OGR

  • Pakiety rgdal, raster czy sf wykorzystują biblioteki GDAL/OGR w R do wczytywania i zapisywania danych przestrzennych
  • http://www.gdal.org/
  • GDAL to biblioteka zawierająca funkcje służące do odczytywania i zapisywania danych w formatach rastrowych
  • OGR to biblioteka służąca to odczytywania i zapisywania danych w formatach wektorowych

2.1.5 PROJ

  • Pakiety sp, raster, czy sf używają biblioteki PROJ do określania i konwersji układów współrzędnych
  • https://proj4.org/
  • Dane przestrzenne powinny być zawsze powiązane z układem współrzędnych
  • PROJ - to biblioteka pozwalająca na identyfikację oraz konwersję pomiędzy różnymi układami współrzędnych
  • Strona http://www.spatialreference.org/ zawiera bazę danych układów współrzędnych
  • Układy współrzędnych mogą być opisywane na szereg sposobów, np poprzez:

    • Reprezentację proj4string - określa ona szereg własności ukladu współrzędnych
    • Kod EPSG (ang. European Petroleum Survey Group) - pozwala on na łatwe identyfikowanie układów współrzędnych.
  • Przykładowo, układ PL 1992 może być określony jako:

"+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"

  • …lub też za pomocą kodu EPSG:

"+init=epsg:2180"

2.1.6 Układ geograficzny

  • Proporcje pomiędzy współrzędną oznaczającą długość geograficzną (X) a oznaczającą szerokość geograficzną (Y) nie są równe 1:1
  • Wielkość oczka siatki jest zmienna
  • Nie pozwala to na proste określanie odległości czy powierzchni
  • Jednostka mapy jest abstrakcyjna
  • Powyższe cechy układów geograficznych powodują, że do większości algorytmów w geostatystyce wykorzystywane są układy współrzędnych prostokątnych płaskich

2.1.7 Układ współrzędnych prostokątnych płaskich

  • Określane są w miarach liniowych (np. metrach)
  • Proporcje między współrzędną X a Y są równe 1:1
  • Wielkość oczka jest stała

2.1.8 GEOS

  • Pakiety rgeos czy sf używają biblioteki GEOS do wykonywania operacji przestrzennych
  • http://geos.osgeo.org/
  • GEOS to biblioteka pozwalająca na przetwarzanie obiektów przestrzennych
  • Przykładowe funkcje tej biblioteki to tworzenie buforów, wyliczanie centroidów, określanie relacji topologicznych (np. przecina, zawiera, etc.) i wiele innych

2.2 Import danych

R pozwala na odczytywanie danych przestrzennych z wielu formatów. Do najpopularniejszych należą dane tekstowe z plików .csv, dane wektorowe z plików .shp, dane rastrowe z plików w formacie GeoTIFF, oraz bazy danych przestrzennych z plików .gpkg.

2.2.1 Format .csv (dane punktowe)

Dane z plików tekstowych (np. .csv) można odczytać za pomocą uogólnionej funkcji read.table() lub też funkcji szczegółowych - read.csv() lub read.csv2().

##       srtm clc      temp      ndvi      savi        x        y
## 1 175.7430   1 13.852222 0.6158061 0.4189449 750298.0 716731.6
## 2 149.8111   1 15.484209 0.5558816 0.3794864 753482.9 717331.4
## 3 272.8583  NA 12.760814 0.6067462 0.3745572 747242.5 720589.0
## 4 187.2777   1 14.324648 0.3756170 0.2386246 755798.9 718828.1
## 5 260.1366   1 15.908549 0.4598393 0.3087599 746963.5 717533.5
## 6 160.1416   2  9.941118 0.5600288 0.3453627 756801.6 720474.1

Po wczytaniu za pomocą funkcji read.csv(), nowy obiekt (np. punkty_sp) jest reprezentowany za pomocą klasy nieprzestrzennej data.frame. Aby obiekt został przetworzony do klasy przestrzennej, konieczne jest nadanie mu współrzędnych. W tym wypadku współrzędne znajdowały się w kolumnach x oraz y. Nadanie układu współrzędnych odbywa się poprzez funkcję coordinates().

## Object of class SpatialPointsDataFrame
## Coordinates:
##        min      max
## x 745592.5 756967.8
## y 712642.4 721238.7
## Is projected: NA 
## proj4string : [NA]
## Number of points: 248
## Data attributes:
##       srtm            clc             temp             ndvi       
##  Min.   :146.5   Min.   :1.000   Min.   : 7.883   Min.   :0.2024  
##  1st Qu.:191.5   1st Qu.:1.000   1st Qu.:12.003   1st Qu.:0.4636  
##  Median :217.9   Median :1.000   Median :14.941   Median :0.5154  
##  Mean   :214.9   Mean   :1.481   Mean   :15.273   Mean   :0.5047  
##  3rd Qu.:239.5   3rd Qu.:2.000   3rd Qu.:17.630   3rd Qu.:0.5742  
##  Max.   :278.4   Max.   :4.000   Max.   :24.945   Max.   :0.6597  
##  NA's   :3       NA's   :5       NA's   :1        NA's   :1       
##       savi       
##  Min.   :0.0824  
##  1st Qu.:0.2935  
##  Median :0.3256  
##  Mean   :0.3176  
##  3rd Qu.:0.3594  
##  Max.   :0.4404  
##  NA's   :1

Ważne, ale nie wymagane, jest także dodanie informacji o układzie przestrzennym danych za pomocą funkcji proj4string().

## Object of class SpatialPointsDataFrame
## Coordinates:
##        min      max
## x 745592.5 756967.8
## y 712642.4 721238.7
## Is projected: TRUE 
## proj4string :
## [+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 +no_defs]
## Number of points: 248
## Data attributes:
##       srtm            clc             temp             ndvi       
##  Min.   :146.5   Min.   :1.000   Min.   : 7.883   Min.   :0.2024  
##  1st Qu.:191.5   1st Qu.:1.000   1st Qu.:12.003   1st Qu.:0.4636  
##  Median :217.9   Median :1.000   Median :14.941   Median :0.5154  
##  Mean   :214.9   Mean   :1.481   Mean   :15.273   Mean   :0.5047  
##  3rd Qu.:239.5   3rd Qu.:2.000   3rd Qu.:17.630   3rd Qu.:0.5742  
##  Max.   :278.4   Max.   :4.000   Max.   :24.945   Max.   :0.6597  
##  NA's   :3       NA's   :5       NA's   :1        NA's   :1       
##       savi       
##  Min.   :0.0824  
##  1st Qu.:0.2935  
##  Median :0.3256  
##  Mean   :0.3176  
##  3rd Qu.:0.3594  
##  Max.   :0.4404  
##  NA's   :1

Proste wyświetlenie uzyskanych danych klasy przestrzennej, np. SpatialPointsDataFrame, można uzyskać za pomocą funkcji plot().

2.2.2 Dane poligonowe (formaty gisowe)

Dane wektorowe (np. shapefile) można odczytać za pomocą funkcji readOGR() z pakietu rgdal. Dla danych w formacie shapefile, przyjmuje ona co najmniej dwa argumenty - dsn oraz layer. Argument dsn określa folder, w którym znajdują się dane. W przypadku, gdy dane znajdują się w folderze roboczym należy ten argument określić za pomocą znaku kropki ("."). Argument layer wymaga podania nazwy pliku bez jego rozszerzenia. Przykładowo, gdy pliki nazywają się granica.dbf, granica.prj, granica.shp, oraz granica.shx - konieczne jest podanie w argumencie layer jedynie nazwy granica.

## OGR data source with driver: GPKG 
## Source: "/home/travis/build/Nowosad/geostat_book/dane/granica.gpkg", layer: "granica.gpkg"
## with 1 features
## It has 3 fields

2.2.3 Rastry

Istnieje kilka sposobów odczytu danych rastrowych w R. Do najpopularniejszych należą funkcje readGDAL() z pakietu rgdal oraz raster() z pakietu raster. Należy w nich jedynie podać ścieżkę do pliku rastrowego.

## dane/siatka.tif has GDAL driver GTiff 
## and has 96 rows and 127 columns

2.3 Przeglądanie danych przestrzennych

2.3.1 Struktura obiektu

Podstawowe informacje o obiektach można uzyskać za pomocą funkcji str():

## Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
##   ..@ data       :'data.frame':  248 obs. of  5 variables:
##   .. ..$ srtm: num [1:248] 176 150 273 187 260 ...
##   .. ..$ clc : int [1:248] 1 1 NA 1 1 2 2 1 1 2 ...
##   .. ..$ temp: num [1:248] 13.9 15.5 12.8 14.3 15.9 ...
##   .. ..$ ndvi: num [1:248] 0.616 0.556 0.607 0.376 0.46 ...
##   .. ..$ savi: num [1:248] 0.419 0.379 0.375 0.239 0.309 ...
##   ..@ coords.nrs : int [1:2] 6 7
##   ..@ coords     : num [1:248, 1:2] 750298 753483 747243 755799 746964 ...
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:248] "1" "2" "3" "4" ...
##   .. .. ..$ : chr [1:2] "x" "y"
##   ..@ bbox       : num [1:2, 1:2] 745593 712642 756968 721239
##   .. ..- 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,"| __truncated__

Obiekty klas Spatial* mają pięć elementów (ang. slots) rozpoczynających się od symbolu @:

  • @data - tabela atrybutów
  • @coords.nrs - numer kolumn zawierających współrzędne
  • @coords - współrzędne kolejnych elementów
  • @bbox - obwiednia
  • @proj4string - definicja układu współrzędnych

2.3.2 Tabla atrybutów

@data jest obiektem klasy data.frame zawierającym informacje o kolejnych punktach w tym zbiorze.

##       srtm clc      temp      ndvi      savi
## 1 175.7430   1 13.852222 0.6158061 0.4189449
## 2 149.8111   1 15.484209 0.5558816 0.3794864
## 3 272.8583  NA 12.760814 0.6067462 0.3745572
## 4 187.2777   1 14.324648 0.3756170 0.2386246
## 5 260.1366   1 15.908549 0.4598393 0.3087599
## 6 160.1416   2  9.941118 0.5600288 0.3453627

2.3.3 Współrzędne

Element @coords.nrs określa numer kolumn zawierających współrzędne w oryginalnym zbiorze danych.

## [1] 6 7

Kolejny element, @coords definiuje współrzędne kolejnych elementów zadanego obiektu.

##          x        y
## 1 750298.0 716731.6
## 2 753482.9 717331.4
## 3 747242.5 720589.0
## 4 755798.9 718828.1
## 5 746963.5 717533.5
## 6 756801.6 720474.1

2.3.4 Obwiednia

Element @bbox określa zasięg przestrzenny danych w jednostkach mapy.

##        min      max
## x 745592.5 756967.8
## y 712642.4 721238.7

2.3.5 Układ współrzędnych

@proj4string reprezentuje definicję układu współrzędnych.

## CRS arguments:
##  +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 +no_defs

2.4 Eksport danych

2.4.1 Zapisywanie danych wektorowych

R pozwala również na zapisywanie danych przestrzennych. W przypadku zapisu danych wektorowych za pomocą funkcji writeOGR() konieczne jest podanie nazwy zapisywanego obiektu (np. granica), folderu w którym chcemy zapisać plik (np. dane), nazwę zapisywanych plików bez rozszerzenia (np. nowa_granica), oraz sterownik - w przypadku danych shapefile jest to ESRI Shapefile.

Zapisywanie w bazach danych przestrzennych odbywa się poprzez nazwy zapisywanego obiektu (np. granica), ścieżki do zapisywanego pliku (np. dane/nowa_granica.gpkg), nazwę dla nowoutworzonej warstwy (np. granica), oraz sterownik - w przypadku danych GeoPackage jest to GPKG.

2.4.2 Zapisywanie danych rastrowych

Funkcja writeGDAL() pozwala na zapisywanie danych rastrowych. Wymaga ona podania dwóch argumentów - nazwy zapisywanego obiektu (np. siatka_sp) oraz ścieżki i nazwy nowego pliku wraz z rozszerzeniem (np. dane/nowy_siatka.tif).

Możliwe jest także użycie szeregu dodatkowych argumentów, np.: - type - określającego typ danych wyjściowych - options - pozwalającego na użycie dodatkowych opcji sterownika, np. użycie kompresji danych

2.5 Wizualizacja danych 2D

Do wizualizacji danych przestrzennych w R służy co najmniej kilkanaście różnych pakietów. Poniżej pokazane są przykłady kilku najprostszych funkcji - plot() oraz spplot() z pakietu sp oraz levelplot() z pakietu rasterVis. Więcej o wizualizacji danych przestrzennych można przeczytać w rozdziale Making maps with R książki Geocomputation with R.

2.5.1 Dane punktowe

Funkcja plot() idealnie nadaje się do szybkiego przyjrzenia się, np. rodzajowi próbkowania danych.

Funkcja spplot() w prosty sposób pozwala na obejrzenie rozkładu wartości wybranej zmiennej. Należy w niej podać nazwę obiektu oraz nazwę wyświetlanej zmiennej. Poniżej można zobaczyć przykłady dla zmiennych temp oraz srtm.

2.5.2 Dane punktowe - kategorie

Nie zawsze dane mają postać ciągłych wartości - bywają one również określeniami różnych klas. W takich sytuacjach należy wcześniej przetworzyć typ danych do postaci kategorycznej (as.factor()). Następnie można je wyświetlić za pomocą funkcji spplot().

2.5.3 Rastry

Wyświetlanie danych w formacie rastrowym może odbyć się z użyciem funkcji plot() lub spplot(). Pierwsza z nich wymaga zdefiniowania obiektu do wizualizacji i domyślnie wyświetla pierwszą warstwę obiektu. Druga z nich domyślnie wyświetla wszystkie dostępne warstwy lub też pojedynczą warstwę poprzez podanie jej nazwy:

2.6 Zadania

  1. Wczytaj dane z pliku punkty_ndvi.csv i przetwórz je do postaci obiektu przestrzennego ndvi_p.
  2. Stwórz mapę zmiennej ndvi z nowoutworzonego obiektu. (Dodatkowo: zapisz mapę do pliku graficznego punkty_ndvi.png).
  3. Zapisz obiekt ndvi_p do formatu GeoPackage oraz do formatu Esri Shapefile.