Wstęp

Tworzenie rozbudowanych narzędzi, których zadaniem jest przetwarzanie zbiorów danych geo-przestrzennych wiąże się z koniecznością przewidzenia i przeciwdziałania wielu sytuacjom, które mogą mieć miejsce, (a nawet na pewno wystąpią!) w trakcie przetwarzania danych wejściowych. Działanie skryptu, ze względu na objętość danych geo-przestrzennych jest wolne (może trwać godziny, a nawet dni), dlatego wszelkie problemy powinny być obsłużone zanim rozpocznie się właściwa część skryptu. Pół biedy, jeżeli w wyniku problemów skrypt się nie uruchomi. Można wtedy usunąć przyczyny i uruchomić go ponownie. Duże gorzej jeżeli się uruchomi i zwróci błędne wyniki, lub brak wyników. Należy tu zaznaczyć, że omawiane problemy nie dotyczą błędów w algorytmach, ale problemów które przychodzą wraz z danymi lub użytkownikiem. Część tych problemów leży po stronie użytkownika skryptu ale część ma charakter obiektywny i wynika z charakteru danych. Przyczyny problemów mogą być następujące:

Nieprawidłowy typ danych
Niepoprawne dane wejściowe i/lub parametry
Niezgodność georeferencji
Niezgodność przestrzenna
Braki w danych
Wyjątki (naturalne sytuacje, które nie mogą być obsłużone przez algorytm)

R posiada rozbudowane biblioteki analizy danych geo-przestrzennych, które powinny (a nie zawsze to robią) zadbać o prawidłowe obsłużenie tych sytuacji. Im bardziej rozbudowany skrypt, ale też im bardziej ogólne narzędzia stosujemy, tym częściej obsługa powyższych sytuacji spada na nas – autorów. Najprostszą sytuacją jest zatrzymanie skryptu, i poinformowanie użytkownika o problemach, ale dotyczy to jedynie najprostszych sytuacji. Często zachodzi sytuacja, że brakuje danej i następuje próba wykonania operacji na NA lub operacji niedozwolonej np. log(-1). W takiej sytuacji powinniśmy odpowiednio zareagować, aby z jednej strony nie przerwać obliczeń a z drugiej zaznaczyć taki fakt w wynikach.

Wymagane biblioteki: raster, sp, geos, argparser
Wymagana wiedza: podstawy programowania

Spójność danych wejściowych

Ta część oceny na na celu sprawdzenie, czy dane mogą być w ogóle przetwarzane. Obejmuje to sprawdzenie typu danych, zakresu i wielkości, georeferencji i zasięgu przestrzennego.

Typy danych

Typ danych sprawdzamy funkcją class. Możemy porównać zwróconą wartość (tekst) z nazwą wymaganej klasy i jeżeli nie będzie spełniona, to zgłaszamy błąd i zatrzymujemy funkcję poleceniem stop().

library(raster)

## Loading required package: sp

library(rgdal)

## 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

rast <- raster("DATA/pop_low.tif")
vect <- readOGR(dsn="DATA/simple",layer="points")

## OGR data source with driver: ESRI Shapefile 
## Source: "DATA/simple", layer: "points"
## with 9 features
## It has 2 fields

geoprocessing <- function (x) {
  if(class(x) != "SpatialPointsDataFrame") stop(paste("wymagany: SpatialPointDataFrame, a jest:",class(x)))
  message("Udało się!")
}
geoprocessing(3)

## Error in geoprocessing(3): wymagany: SpatialPointDataFrame, a jest: numeric

geoprocessing(rast)

## Error in geoprocessing(rast): wymagany: SpatialPointDataFrame, a jest: RasterLayer

geoprocessing(vect)

## Udało się!

Wczytywanie danych zajmuje dużo czasu (i miejsca w pamięci). Dlatego funkcja readOGR() pozwala sprawdzić typ geometrii zanim zacznie wczytywać dane (nie robi tego dobrze ani przejrzyście). Dlatego jeszcze bardziej sensowym rozwiązaniem jest sprawdzenie metadanych poznanymi wcześniej funkcjami: ogrinfo() i gdalInfo().

vect <- readOGR(dsn="DATA/simple",layer="lines",require_geomType = "wkbPoint") #sprawdzenie czy punkty. zgłasza błąd

## Error in ogrInfo(dsn = dsn, layer = layer, encoding = encoding, use_iconv = use_iconv, : wkbPointnot inwkbLineString

read_points <- function(dsn,layer) {  ##TO JEST ZLE##
  vect_meta <- ogrInfo(dsn=dsn,layer=layer) #pobieramy metadane
  if(vect_meta$eType != 2) stop(paste("wymagane: points, a jest:",vect_meta$eType))
  readOGR(dsn=dsn,layer=layer)
}

vect <- read_points(dsn="DATA/simple",layer="points")

## Error in read_points(dsn = "DATA/simple", layer = "points"): wymagane: points, a jest: 1

vect <- read_points(dsn="DATA/simple",layer="lines")

## OGR data source with driver: ESRI Shapefile 
## Source: "DATA/simple", layer: "lines"
## with 1 features
## It has 1 fields

Sprawdzanie poprawności danych wejściowych

Sprawdzenie poprawności danych wejściowych jest następnym krokiem. Metadane zawracane przez funkcję ogrInfo() i GDALinfo() pozwalają sprawdzić wiele aspektów danych zanim rozpoczniemy je wczytywać. Metadane dla danych wektorowych:

vect_meta <- ogrInfo(dsn="DATA/simple",layer="points") 
str(vect_meta)

## List of 12
##  $ nrows        : int 9
##  $ nitems       : int 2
##  $ iteminfo     :List of 5
##   ..$ name        : chr [1:2] "id" "id2"
##   ..$ type        : int [1:2] 12 12
##   ..$ length      : int [1:2] 10 10
##   ..$ typeName    : chr [1:2] "Integer64" "Integer64"
##   ..$ maxListCount: int [1:2] 0 0
##  $ driver       : chr "ESRI Shapefile"
##  $ extent       : num [1:4] 956936 1846136 996117 1863321
##  $ nListFields  : int 0
##  $ have_features: logi TRUE
##  $ eType        : int 1
##  $ with_z       : int 0
##  $ dsn          : chr "DATA/simple"
##  $ layer        : chr "points"
##  $ p4s          : 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 "
##  - attr(*, "LDID")= int 0
##  - attr(*, "class")= chr "ogrinfo"

najważniejsze metadane obiektów wektorowych to:

nrows - liczba wierszy w tabeli atrybutów
nitems - liczba kolumn w tabeli
iteminfo - informacje o tabeli (typy danych, nazwy kolumn, brak informacji o wartościach)
extent - zasięg przestrzenny
layer - typ danych
p4s - projekcja
driver - typ danych
have_features - ważny parametr, pozwala unikać pustych plików

Niestety, metadane rastrowe są bardziej skompilowane, a dostęp do nich odbywa się poprzez nazwy pól… albo poprzez atrybuty:

rast_meta <- GDALinfo("DATA/pop_low.tif") 
str(rast_meta)

## Class 'GDALobj'  atomic [1:9] 127 201 1 953085 1838625 ...
##   ..- attr(*, "ysign")= num -1
##   ..- attr(*, "driver")= chr "GTiff"
##   ..- attr(*, "projection")= 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 "
##   ..- attr(*, "file")= chr "DATA/pop_low.tif"
##   ..- attr(*, "df")='data.frame':    1 obs. of  9 variables:
##   .. ..$ GDType        : Factor w/ 1 level "Float64": 1
##   .. ..$ hasNoDataValue: logi TRUE
##   .. ..$ NoDataValue   : num -1.7e+308
##   .. ..$ blockSize1    : int 5
##   .. ..$ blockSize2    : int 201
##   .. ..$ Bmin          : num 0
##   .. ..$ Bmax          : num 1508
##   .. ..$ Bmean         : num 42.4
##   .. ..$ Bsd           : num 61.6
##   ..- attr(*, "sdf")= logi TRUE
##   ..- attr(*, "mdata")= chr "AREA_OR_POINT=Area"
##   ..- attr(*, "ScaleOffset")= num [1, 1:2] 1 0
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : NULL
##   .. .. ..$ : chr [1:2] "scale" "offset"

names(rast_meta) #nazwy pól dotyczące geometrii

## [1] "rows"      "columns"   "bands"     "ll.x"      "ll.y"      "res.x"    
## [7] "res.y"     "oblique.x" "oblique.y"

statistics <- attr(rast_meta,"df") #data.frame jako atrybut.
str(statistics) #mamy tylko 1 warstwę (band), gdyby było ich więcej to było by więcej wierszy

## 'data.frame':    1 obs. of  9 variables:
##  $ GDType        : Factor w/ 1 level "Float64": 1
##  $ hasNoDataValue: logi TRUE
##  $ NoDataValue   : num -1.7e+308
##  $ blockSize1    : int 5
##  $ blockSize2    : int 201
##  $ Bmin          : num 0
##  $ Bmax          : num 1508
##  $ Bmean         : num 42.4
##  $ Bsd           : num 61.6

attr(rast_meta,"driver") #driver, geoTIFF

## [1] "GTiff"

Geometria - dostępne poprzez nazwy: rows, columns, bands, ll.x, lly, res.x, res.y
Statystyki danych: typ (zgodny ze standardami GDAL), nodata, rozmiar wewnętrzny bloków, min, max, mean, sd

Co można sprawdzać w zakresie poprawności danych? Przede wszystkim należy ograniczać liczbę testów do niezbędnych.

Typ - przede wszystkim jeżeli potrzebujemy dane kategoryzowane, czy typ jest integer lub byte
Czy w danych występują wartości NoData
Jaki jest rozmiar danych (liczba wierszy i kolumn)
Jaka st rozdzielczość danych (czy nie jest zbyt mała lub zbyt duża)
Czy występują wartości ujemne lub zera (Bmin)
Czy w danych występują Inf i -Inf

Sprawdzanie atrybutów wejściowych

W skryptach/funkcjach geo-przetwarzania parametry procesu są przekazywane przez użytkownika skryptu. W przypadku narzędzi udostępnionych, nieodpowiednie parametry mogą doprowadzić do bardzo niekorzystnych skutków, na przykład zbyt dużej liczby iteracji, zużycia zasobów komputera (w połączeniu z rozmiarem danych wejściowych), lub też nie będą mogły zostać przekazane do wbudowanych funkcji (np. ujemny rozmiar). Jak zapobiegać błędom:

O ile to możliwe - należy ustawić wartości parametrów domyślnych
Ustawiania wartości domyślnych funkcji wbudowanych nie należy wymuszać, należy dać swobodę ich ustawienia poprzez …
Pozostałe parametry testować czy ich wartości nie mogą spowodować problemów

Parametry testujemy poprzez funkcję procedurę if oraz funkcję stop(). Poniższy przykład zawiera funkcję z dwoma parametrami swobodnymi, które mają pewne ograniczenia:

parameters <- function(data,niter=10,scale=1) {
  if(niter<2) stop("minimalna liczba powtórzeń to 2 (zalecane 6)")
  if(scale<=0| scale>1) stop("Parametr skali musi być w przedziale (0,1]")
  message("parametry OK, liczę...")
}

parameters(1,niter=6,scale=3) # typ danych nie jest teraz testowany może być cokoliwiek

## Error in parameters(1, niter = 6, scale = 3): Parametr skali musi być w przedziale (0,1]

parameters(1,niter=-1)

## Error in parameters(1, niter = -1): minimalna liczba powtórzeń to 2 (zalecane 6)

parameters(1) #wartości domyślne

## parametry OK, liczę...

Praca z georeferencjami

W przypadku gdy dane przekazane do skryptu mają różną projekcję, może dojść do sytuacji, że obliczenia zostaną wykonane ale wyniki obliczeń będą (rażąco) błędne. Jednakże ze względu na niedoskonałości współpracy pomiędzy różnymi aplikacjami a GDAL/OGR pliki o tej samej projekcji, mogą mieć nieznacznie różne reprezentacje (sposób zapisu) PROJ4 i sprawiać wrażenie różnych, pomimo że nie będą.

require(raster) 
require(rgdal)

proc_ref <- CRS("+init=epsg:5070") # docelowa projekcja, ta sama którą posiadają wczytywane dane
raster_ref <- GDALSpatialRef("DATA/pop_low.tif")
proc_ref

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

raster_ref #różnią się drobiazgami ale główne elementy są takie same

## [1] "+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 "

compareCRS(proc_ref,raster_ref) # zwłykłe porownanie, pokazuje że to ta sama referencja

## [1] TRUE

compareCRS(proc_ref,raster_ref,verbatim = TRUE) # pełne, że nie, literalne sprawdzenie

## [1] FALSE

#======
layer <- brick("DATA/pop_low.tif") 
layer_ref <- crs(layer)
layer_ref

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

compareCRS(proc_ref,layer_ref)

## [1] TRUE

compareCRS(raster_ref,layer_ref) #wszystkie są jednakowe.

## [1] TRUE

#======
proc_ref <- CRS("+init=epsg:4326") # docelowa projekcja, inna niż ta którą posiadają wczytywane dane
compareCRS(proc_ref,layer_ref)

## [1] FALSE

compareCRS(proc_ref,layer_ref,verbatim = TRUE) # oba fałszywe

## [1] FALSE

Testy zgodności geo-referencji danych rastrowych należy przeprowadzić przed ich wczytywaniem. Negatywny wynik testu nie oznacza jednak konieczności zakończenia pracy, a jedynie konieczność ich reprojekcji. Należy jednak pamiętać że reprojekcja dużego zbioru danych rastrowych mogą być czasochłonna i zużyć zasoby pamięci. W takim przypadku należy przerwać działanie skryptu i nakazać ujednolicenie projekcji przed ich przekazaniem do skryptu. Pozostaje problem, który ze zbiorów danych ma właściwą referencję? W praktyce powinien to być główny zbiór danych (przekazywany do funkcji jako pierwszy parametr). Dodatkowo, jeżeli te same projekcje zapisane są przy pomocy różnych łańcuchów WKT, należy ujednolicić je do jednego formatu.

rasterfile1 <- ("DATA/pop_low.tif")
rasterfile2 <- ("DATA/landcovWGS.tif")

process_geodata <- function(rasterfile1,rasterfile2,force.reprojection=FALSE) { # +możliwość wymuszania reprojekcji
  
  reproject=FALSE # nie robimy reprojekcji
  if(compareCRS(GDALSpatialRef(rasterfile1),GDALSpatialRef(rasterfile2))==FALSE) {
    if(force.reprojection==TRUE) {
      reproject <- TRUE # użytkownik sobie życzy reprojekcji
    } else {
      stop("Różne układy odniesienia w danych wejściowych") #test zgodności ukladów odniesienia 
    }
  }
  
  raster1 <- raster(rasterfile1) #dopiero teraz wczytujemy dane
  raster2 <- raster(rasterfile2)
  
  if(reproject) { 
    message("wymuszono reprojekcję, proszę czekać....")
    raster2 <- projectRaster(raster2,raster1) # używamy projekcji raster1 jako docelowej 
  }
  
  
  if(!compareCRS(raster1,raster2,verbatim = TRUE)) { #sprawdzamy czy oba crs-y mają dokładnie ten sam kod
    message(paste("kopiuję referencję z",rasterfile1,"do",rasterfile2))
    crs(raster2) <- crs(raster1) # ujednolicamy kody projekcji, to nie jest reprojekcja
  }
  
  message("gotowe, liczę...") #dalsze obliczenia
}  

process_geodata("DATA/pop_low.tif","DATA/landcovWGS.tif") # różne projekcje, brak wymuszenia reprojekcji

## Error in process_geodata("DATA/pop_low.tif", "DATA/landcovWGS.tif"): Różne układy odniesienia w danych wejściowych

process_geodata("DATA/pop_low.tif","DATA/landcovWGS.tif",force.reprojection = TRUE) # różne prokekcje, wymuszenie reprojekcji

## wymuszono reprojekcję, proszę czekać....

## gotowe, liczę...

process_geodata("DATA/pop_low.tif","DATA/landcov_low.tif") # jednakowe projekcje

## gotowe, liczę...

W przypadku reprojekcji danych wektorowych, sposób postępowania jest bardzo podobny, z tą różnicą, że jeżeli używamy równocześnie danych rastrowych i wektorowych ze względu na różnice w wydajności i zniekształcenia, projekcji powinny podlegać dane wektorowe, ewentualnie można później przywrócić oryginalną projekcję.

require(raster)
require(rgdal)

rast <- raster("DATA/pop_low.tif")
vect <- readOGR("DATA/simple","lines")

## OGR data source with driver: ESRI Shapefile 
## Source: "DATA/simple", layer: "lines"
## with 1 features
## It has 1 fields

crs(rast)

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

crs(vect) # różne na poziomie verbatim

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

process_geodata <- function(rast,vect) {
  
  orginal_vect_reference <- crs(vect)
  if(!compareCRS(crs(rast),crs(vect))) {
    message("Wykonuję reprojekcę")
    vect <- spTransform(vect,crs(rast))
  }
  

  if(!compareCRS(crs(rast),crs(vect),verbatim = TRUE)) {
     message("kopiuję referencję z rast do vect")
    crs(vect) <- crs(rast) # skopiowanie projekcji
  }
  
  message("obliczenia...")
  result_vect <- vect #rezultat jakiś obliczeń
  if(!compareCRS(crs(result_vect),orginal_vect_reference)) {
    message("przywracam oryginalną referencję")  
    crs(result_vect) <- orginal_vect_reference  
  }
  result_vect
}

vect2 <- process_geodata(rast,vect)

## kopiuję referencję z rast do vect

## obliczenia...

Analiza zgodności geoprzestrzennej

Analiza zgodności geo-przestrzennej to kolejny krok jaki należy wykonać przy przetwarzaniu więcej niż jednego zbioru danych geo-przestrzennych. Aby wykonać chociażby popularne operacje algebry map warstwy muszą się idealnie nakładać. W przypadku danych wektorowych sytuacja jest często jeszcze trudniejsza - dane nie muszą się nakładać (i z reguły nie będą), ale z drugiej strony powinny należeć do jednego obszaru geograficznego, aby analiza miała sens.

Do sprawdzenia zgodności danych rastrowych służy funkcja compareRast(), która pozwala na porównanie dwóch zbiorów rastrowych. Domyślnie funkcja dokonuje pełnego porównania obejmującego rozdzielczość, zasięg (extent), crs, rozdzielczość, liczbę wierszy i kolumn. Opcjonalnie funkcja może sprawdzić również wartości. Uwaga!!!. Funkcja działa w sposób podobny do stopifnot() to znaczy zatrzymuje działanie programu, zgłaszając błąd. Jeżeli chcemy kontynuować pracę (na przykład wykonując dostosowanie danych) należy ustawić parametr stopiffalse na FALSE (domyślnie jest TRUE).

library(raster)
rast_pop <- raster("DATA/pop_low.tif")
rast_land <- raster("DATA/landcov_low.tif")
rast_WGS <- raster("DATA/landcovWGS.tif")
compareRaster(rast_pop,rast_land)

## Error in compareRaster(rast_pop, rast_land): different extent

compareRaster(rast_pop,rast_land,stopiffalse = FALSE) #wiemy że nie

## [1] FALSE

rast_pop_agg <- disaggregate(rast_pop,fact=2) #zwiększamy rozdzielczość
compareRaster(rast_pop,rast_pop_agg,rowcol=FALSE) # raster o takim samym zasięgu ale innej rozdzielczosci

## [1] TRUE

W przypadku, gdy dwa pliki rastrowe mają różną geometrię (na przykład extent, rozdzielczość) możemy zmodyfikować nowy raster do wymaganych parametrów. Służy do tego funkcje resample() i projectRaster(). ProjectRaster zintegruje geometrię rastrów tylko w sytuacji, gdy jako referencję podamy obiekt raster a nie crs. Należy zauważyć, że w zależności od wybranej metody przetwarzania danych otrzymamy nieco inne dane wyjściowe. Wynika to z faktu, że mamy inną kolejność postępowania: najpierw reprojekcja potem resampling lub odwrotnie.

rast_land_rep <- projectRaster(rast_WGS,rast_pop,method = "ngb") #reprojekcja 4326 do 5070
compareRaster(rast_pop,rast_land_rep)

## [1] TRUE

rast_land_res <- resample(rast_land,rast_pop,method = "ngb")
compareRaster(rast_land_res,rast_land_rep) #geometria się zgadza się

## [1] TRUE

compareRaster(rast_land_res,rast_land_rep,values=TRUE) #nie zgadzają się niektóre wartości

## Error in compareRaster(rast_land_res, rast_land_rep, values = TRUE): not all objects have the same values

rast_land_rep2 <- projectRaster(rast_land,rast_pop,method = "ngb")
compareRaster(rast_land_res,rast_land_rep2) #geometria się zgadza się

## [1] TRUE

compareRaster(rast_land_res,rast_land_rep2,values=TRUE) #zgadzają się wartości

## [1] TRUE

Porównywanie zakresu/typu obiektów wektorowych [i rastrowych] ma zupełnie inne zadanie: przede wszystkim pozwala unikać sytuacji przetwarzania danych, które nie mogą dać pozytywnego wyniku (na przykład próbkowanie, gdy zasięg pliku wektorowego nie pokrywa się z plikiem rastrowym). Taki test możemy przeprowadzić wykorzystując funkcje z biblioteki rgeos: gWithin() testujących czy jeden bbox zawiera się w drugim, czy gOverlaps(), która sprawdza czy daw bboxy się częściowo pokrywają. Bounding box musi być wcześniej zamieniony w poligon, co można wykonać poleceniami extent() i as.Polygon() z biblioteki raster. Aby nie stosować dwóch funkcji do testowania czy dwa obiekty wektorowe współdzielą ze sobą zasięg można zastosować zaprzeczenie funkcji gDisjoint(), która zwraca prawdę, w przypadku gdy obiekty nie mają ze sobą nic wspólnego. W przypadku odległych obiektów można zastosować funkcję gDistance(), ale jej wynik trzeba porównać z pewną założoną minimalną wartością progową.

require(rgdal)
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: TRUE

rast_pop <- raster("DATA/pop_low.tif")
vect_lines <- readOGR("DATA/simple","lines")

## OGR data source with driver: ESRI Shapefile 
## Source: "DATA/simple", layer: "lines"
## with 1 features
## It has 1 fields

vect_points <- readOGR("DATA/simple","points")

## OGR data source with driver: ESRI Shapefile 
## Source: "DATA/simple", layer: "points"
## with 9 features
## It has 2 fields

pol_raster <- as(extent(rast_pop),"SpatialPolygons")
pol_vector_lines <- as(extent(bbox(vect_lines)),"SpatialPolygons")
pol_vector_points <- as(extent(bbox(vect_points)),"SpatialPolygons")
plot(gDifference(pol_raster,pol_vector_lines),col=2) #zasięgi bboxów: rastra (zewnętrzny) i linii(wewnętrzny)

gWithin(pol_vector_lines,pol_raster) # wektor w rastrze

## [1] TRUE

gOverlaps(pol_vector_lines,pol_raster) #overlaps wymaga częściowego pokrycia

## [1] FALSE

gOverlaps(pol_vector_lines,pol_vector_points) #bounding boxy pokrywają się częsciowo

## [1] TRUE

!gDisjoint(pol_vector_lines,pol_raster) #funkcja, która testuje czy elementy nie mają ze sobą nic wspólnego, z not na początku uogólnia  zarówno częściowe jak i całkowite nakładanie

## [1] TRUE

Pisanie skryptów i funkcji geoprzetwarzania

Dzięki licznym funkcjom geo-przetwarzania oraz dobrej współpracy z zewnętrznymi systemami GIS R bardzo dobrze nadaje się do pisania zarówno funkcji jak i skryptów geo-przetwarzania. Funkcje geo-przetwarzania wymagają uruchomienia środowiska R i wywołania już napisanych funkcji, natomiast skrypty geo-przetwarzania to gotowe procedury, które użytkownik może traktować jak programy uruchamiane z linii poleceń, bez znajomości ich struktury wewnętrznej W tej części lekcji zapoznamy się jak pisać i wywoływać takie procedury.

Analizę przeprowadzimy na podstawie prostej funkcji/ skryptu, która dokonuje próbkowania rastra przy pomocy mapy punktów. W przypadku jej braku, skrypt wygeneruje losowy zbiór danych.

Tworzenie funkcji

Będziemy potrzebowali dwie funkcje: sample_random() i sample_by_points(), funkcja pierwsza jako argument przyjmuje obiekt raster i liczbę punktów, druga jako argument przyjmuje raster i shapefile. Obie jako wynik zwracają SpatialPointDataFrame. Funkcje muszą sprawdzić, czy obiekty się zgadzają, czy mają prawidłowy CRS, oraz czy zasięg obiektów się zgadza. Jeżeli nie powinny przerwać działanie. istotą funkcji w R jest to żeby obsługiwały obiekty R a nie zewnętrzne obiekty geo-przestrzenne. Konwersja powinna być wykonana przed przekazaniem argumentów do funkcji. W przeciwnym wypadku funkcje nie będą uniwersalne, tj. nie będzie można ich wykorzystywać na bardziej ogólnym poziomie. Funkcje powinny wywoływać biblioteki poleceniem require(), które doczyta brakujące biblioteki.

# funkcja losowa
sample_random <- function(rast,number=1) {
  require(raster)
  if(class(rast)!="RasterLayer") stop(paste("wymagany: Raster, a jest:",class(rast)))
  if(number<1) stop ("wymagany minimum 1 punkt")
  sampleRandom(rast,size=number,sp=TRUE)
}
input_raster <- raster("DATA/pop_low.tif")
output <- sample_random(input_raster,10)
str(output)

## Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
##   ..@ data       :'data.frame':  10 obs. of  1 variable:
##   .. ..$ pop_low: num [1:10] 4.42 38.11 0 23.69 29.81 ...
##   ..@ coords.nrs : num(0) 
##   ..@ coords     : num [1:10, 1:2] 997365 981045 988245 986565 997605 ...
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : NULL
##   .. .. ..$ : chr [1:2] "x" "y"
##   ..@ bbox       : num [1:2, 1:2] 957765 1845945 997605 1866345
##   .. ..- 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"

Funkcja z wektorem jest dużo bardziej skomplikowana, ponieważ zawiera kilka testów oraz konwersji.

# funkcja na podstawie pliku
sample_by_points <- function(rast,vect) {
  require(raster)
  require(rgeos)

  if(class(rast)!="RasterLayer") stop(paste("wymagany: Raster, a jest:",class(x)))
  if(class(vect) != "SpatialPointsDataFrame") stop(paste("wymagany: SpatialPointDataFrame, a jest:",class(vect)))
  
  if(!compareCRS(crs(rast),crs(vect))) stop("Obiekty mają różną georeferencję")
  if(!compareCRS(crs(rast),crs(vect),verbatim = TRUE)) {
    warning("kopiuję referencję z rast do vect")
    crs(vect) <- crs(rast) # skopiowanie projekcji
  }
  
  #sprawdzenie zgodności przestrzennej
  pol_raster <- as(extent(rast),"SpatialPolygons")
  pol_vector <- as(extent(extent(bbox(vect))),"SpatialPolygons")
  if(gDisjoint(pol_vector_lines,pol_raster)) stop("Raster i Wektor nie mają części wspólnej")
  if(!gWithin(pol_vector,pol_raster)) warning("Nie wszystkie punkty znajdują się na rastrze")
  
  sp_raster <- as(input_raster,"SpatialGridDataFrame") #wykorzystamy pakiet sp, zamiast funkcji extract
  vect$sample <- over(vect,sp_raster)[[1]]
  vect
}

input_raster <- raster("DATA/pop_low.tif")
input_vector <- shapefile("DATA/simple/points.shp") #alternatywa dla readOGR, ale tylko ESRI SHAPEFILE
output <- sample_by_points(input_raster,input_vector)

## Warning in sample_by_points(input_raster, input_vector): kopiuję referencję
## z rast do vect

Wyniki można zapisać do pliku zewnętrznego lub zachować do dalszego przetwarzania. Natomiast funkcje można zapisać w zewnętrznym pliku (na przykład funkcje.R) i wywoływać poleceniem source().

Parsowanie argumentów wejściowych do skryptu

Jeżeli Użytkownik skryptu ma z niego korzystać bez konieczności analizy jego struktury wewnętrznej. Aby można było obsługiwać argumenty wejściowe w sposób prosty można skorzystać z jednej z bibliotek oferowanych przez R do parsowania argumentów wejściowych do skryptów. W efekcie taki skrypt będzie się zachowywał podobnie jak każdy inny skrypt uruchamiany z linii poleceń.

library(argparser)
p <- arg_parser("skrypt tworzący próbę losową ze zbioru rastrowego", name="program")
p <- add_argument(p,"input",help="raster wejsciowy") #argument pozycyjny, przedostani
p <- add_argument(p,"output",help="wektor wyjściowy")
#argumenty opcjonalne zaznaczone przez --
p <- add_argument(p,"--vector",help="wektor wejsciowy",default=NA) #  default ustawione na NA
p <- add_argument(p,"--column",help="nazwa kolumny",default="sample")
p <- add_argument(p,"--number",help="liczba punktów",type="integer") # argument opcjonalny, integer
p <- add_argument(p,"--overwrite",help="nadpisać",flag=TRUE) #FLAGA
parse_args(p) #zwróci help, bo żadne argumenty nie zostały przekazane

## usage: program [--] [--help] [--overwrite] [--opts OPTS] [--vector VECTOR] [--column COLUMN] [--number NUMBER] input output

## skrypt tworzący próbę losową ze zbioru rastrowego

## positional arguments:

##   input          raster wejsciowy

##   output         wektor wyjściowy

##

## flags:

##   -h, --help         show this help message and exit

##   -o, --overwrite            nadpisać

##

## optional arguments:

##   -x, --opts OPTS            RDS file containing argument values

##   -v, --vector VECTOR            wektor wejsciowy

##   -c, --column COLUMN            nazwa kolumny [default: sample]

##   -n, --number NUMBER            liczba punktów

## Error in parse_args(p): Missing required arguments: expecting 2 values but got ().

Aby polecenia zostały wykonane, należy wywołać funkcję parse arguments, która jako argument pobiera listę argumentów przekazanych do skryptu parse_args(). Argumenty to wektor znakowy (char vector), domyślnie wynik funkcji commandArgs(trailingOnly = TRUE), która zwraca argumenty przekazane do skryptu. Do celów testowych można taką listę przygotować ręcznie. Wektor zawiera następujące po sobie nazwy argumentów i ich wartości. Skrypt zwraca listę zawierającą nazwy argumentów i przypisane wartości. Proszę zwrócić uwagę, że nie podaliśmy argumentu kolumna ale jest on pobrany z wartości default

# ze zględu na konieczność wywołania skryptu w śrdowisku R używam komendy system()
system(" Rscript DATA/parser.R pop_low.tif simple/samples.shp --vector simple/points.shp")

W kolejnych krokach zrezygnujemy z polecenia system i będziemy przekazywać argumenty jako wektor wejściowy, aby obserwować działanie skryptu.

args <- c("DATA/pop_low.tif", # pozycyjny
          "DATA/simple/sample.shp", # pozycyjny
          "--vector","DATA/simple/points.shp",
          "--overwrite") #opcjonalny
parsed <- parse_args(p,args) #to samo co w poprzednim kroku
parsed

## [[1]]
## [1] FALSE
## 
## $help
## [1] FALSE
## 
## $overwrite
## [1] TRUE
## 
## $opts
## [1] NA
## 
## $vector
## [1] "DATA/simple/points.shp"
## 
## $column
## [1] "sample"
## 
## $number
## [1] NA
## 
## $input
## [1] "DATA/pop_low.tif"
## 
## $output
## [1] "DATA/simple/sample.shp"

Parser automatycznie NIE sprawdza czy pliki istnieją, dla niego są to tylko łańcuch znaków. Z tego w pierwszej kolejności należy sprawdzić czy istnieją dane wejściowe i czy dane wyjściowe nie zostaną nadpisane.następnie należy sprawdzić ich zgodność z wymogami (omówione w poprzednich krokach kursu):

#sprawdzenie istenienia plików
if(!file.exists(parsed$input)) stop(paste("nie mogę znaleść pliku", parsed$raster))
if(file.exists(parsed$output) && parsed$overwrite==FALSE) stop(paste("Plik", parsed$raster, "istnieje, użyj --overwrite, aby nadpisać")) #należy zabezpieczyć się przed przypadkowym nadpisaniem

Następny krok to sprawdzenie, który rodzaj opcji wybrał użytkownik: próbkowanie opcjonalnym rastrem lub generowanie własnego zbioru. Test kończy się ustawieniem wartości parsed$number lub zwróceniem komunikatu o błędzie, jeżeli nie podano wektora wejściowego i liczby generowanych elementów. Jeżeli podano oba ZAKŁADAMY że istniejący wektor posiada pierwszeństwo:

# po kolei eliminujemy opcje uniemożliwjaące wykonanie skryptu
if(is.na(parsed$vector) && is.na(parsed$number)) stop("wymagany wektor lub liczba punktów")
if(!is.na(parsed$vector) && !file.exists(parsed$vector)) stop(paste("nie mogę znaleść pliku", parsed$vector)) #jest arguement ale błędna nazwa pliku
if(!is.na(parsed$vector)) 
{
  parsed$number <- NA # wektor jest; number już nas nie obchodzi
} else  {
  if(parsed$number<=0) stop ("liczba punktów number musi wynosić conajmniej 1") #tylko, 
}
# Jeżeli parsed$number ocalał, to znaczy że jego używamy do dalszej pracy

Wczytujemy dane i/lub generujemy losowy plik próbkujący. parsed$number pełni rolę flagi, jest jest NA używamy wektora, jak nie liczby. Funkcje napisane wcześniej wczytujemy poleceniem source().

source("DATA/funkcje.R")
library(raster)
input_raster <- raster(parsed$input)
if(!is.na(parsed$vector)) {
  input_vector <- shapefile(parsed$vector)
  output <- sample_by_points(input_raster,input_vector)
  names(output) <- c(names(input_vector),parsed$column) #nadanie nazwy kolumnie
} else {
  output <- sample_random(input_raster,10)
  names(output) <- parsed$column
}

## Warning in sample_by_points(input_raster, input_vector): kopiuję referencję
## z rast do vect

shapefile(output,file=parsed$output,overwrite=parsed$overwrite) # decyzja o ostatecznym zapisie

Waściwy skrypt geoprzetwarzania

Ze względu na problemy ze ścieżkami w środowisku Rstudio, skrypty należy wywoływać bezpośrednio z linii poleceń

Plik funkcje.R

# funkcja losowa
sample_random <- function(rast,number=1) {
  require(raster)
  if(class(rast)!="RasterLayer") stop(paste("wymagany: Raster, a jest:",class(rast)))
  if(number<1) stop ("wymagany minimum 1 punkt")
  sampleRandom(rast,size=number,sp=TRUE)
}

# funkcja na podstawie pliku
sample_by_points <- function(rast,vect) {
  require(raster)
  require(rgeos)
  
  if(class(rast)!="RasterLayer") stop(paste("wymagany: Raster, a jest:",class(x)))
  if(class(vect) != "SpatialPointsDataFrame") stop(paste("wymagany: SpatialPointDataFrame, a jest:",class(vect)))
  
  if(!compareCRS(crs(rast),crs(vect))) stop("Obiekty mają różną georeferencję")
  if(!compareCRS(crs(rast),crs(vect),verbatim = TRUE)) {
    warning("kopiuję referencję z rast do vect")
    crs(vect) <- crs(rast) # skopiowanie projekcji
  }
  
  #sprawdzenie zgodności przestrzennej
  pol_raster <- as(extent(rast),"SpatialPolygons")
  pol_vector <- as(extent(extent(bbox(vect))),"SpatialPolygons")
  if(gDisjoint(pol_vector_lines,pol_raster)) stop("Raster i Wektor nie mają części wspólnej")
  if(!gWithin(pol_vector,pol_raster)) warning("Nie wszystkie punkty znajdują się na rastrze")
  
  sp_raster <- as(input_raster,"SpatialGridDataFrame") #raster nie posiada funkcji próbkowania na podstawie shp, musimy wykorzystać pakiet sp
  vect$sample <- over(vect,sp_raster)[[1]]
  vect
}

Plik parser.R

# skrypt geoprzetwarzania ======= PARSER ====================================#
library(argparser)
p <- arg_parser("skrypt tworzący próbę losową ze zbioru rastrowego", name="program")
p <- add_argument(p,"input",help="raster wejsciowy") #argument pozycyjny, przedostani
p <- add_argument(p,"output",help="wektor wyjściowy")
#argumenty opcjonalne zaznaczone przez --
p <- add_argument(p,"--vector",help="wektor wejsciowy",default=NA) #  default ustawione na NA
p <- add_argument(p,"--column",help="nazwa kolumny",default="sample")
p <- add_argument(p,"--number",help="liczba punktów",type="integer") # argument opcjonalny, integer
p <- add_argument(p,"--overwrite",help="nadpisać",flag=TRUE) #FLAGA
parsed <- parse_args(p)

#===================================================================================#

#sprawdzenie istenienia plików
if(!file.exists(parsed$input)) stop(paste("nie mogę znaleść pliku", parsed$raster))
if(file.exists(parsed$output) && parsed$overwrite==FALSE) stop(paste("Plik", parsed$raster, "istnieje, użyj --overwrite, aby nadpisać")) #należy zabezpieczyć się przed przypadkowym nadpisaniem

if(is.na(parsed$vector) && is.na(parsed$number)) stop("wymagany wektor lub liczba punktów")
if(!is.na(parsed$vector) && !file.exists(parsed$vector)) stop(paste("nie mogę znaleść pliku", parsed$vector)) 
if(!is.na(parsed$vector)) 
{
  parsed$number <- NA
} else  {
  if(parsed$number<=0) stop ("liczba punktów number musi wynosić conajmniej 1") 
}

#===================================================================================#

source("funkcje.R") # wywołujemy funkcje
library(raster)
input_raster <- raster(parsed$input)
if(!is.na(parsed$vector)) {
  input_vector <- shapefile(parsed$vector)
  output <- sample_by_points(input_raster,input_vector)
  names(output) <- c(names(input_vector),parsed$column) #nadanie nazwy kolumnie
} else {
  output <- sample_random(input_raster,10)
  names(output) <- parsed$column
}

#================================== Zapis wyników ==================================#

shapefile(output,file=parsed$output,overwrite=parsed$overwrite)

Skrypty geo-przetwarzania