Saturday, 25 May 2019

Extraction of raster cell values based on 'rgeos' and 'raster' buffering approaches in R


I am trying to extract cell values from a raster that lie within a certain distance from a given location (i.e. buffering). Evaluating different approaches, I found out that applying gBuffer(..., width = 20000) from package rgeos generally includes fewer cells than extract(..., buffer = 2000) from package raster for one and the same buffer width. Does anybody know what's the difference between these two approaches?


Here's a reproducible example that comes close to my piece of code:



# Required packages
library(dismo)
library(raster)
library(rgeos)

# Sample raster data (and perform reclassification for clarity purposes)
germany <- gmap("Germany")
germany <- reclassify(germany, matrix(c(0, 200, 0,
200, 255, 1), ncol = 3, byrow = TRUE))


# Sample point shape data
nuremberg <- data.frame("lat" = 49.452778, "lon" = 11.077778, "city" = "Nuremberg")
# Set current CRS and reproject to raster CRS
coordinates(nuremberg) <- c("lon", "lat")
projection(nuremberg) <- CRS("+proj=lonlat +ellps=WGS84")
nuremberg <- spTransform(nuremberg, CRS(projection(germany)))

# Extract cell values using extract() only
nbg.val <- unlist(extract(germany, nuremberg, buffer = 20000))


# Extract cell values using gBuffer() and extract()
nbg.bff <- gBuffer(nuremberg, width = 20000)
nbg.val2 <- unlist(extract(germany, nbg.bff))

So far, so good. Now let's compare the results.


# extract() only
table(nbg.val)
nbg.val
0 1
118 91


# gBuffer() and extract()
table(nbg.val2)
nbg.val2
0 1
116 90

Of course, these are minor differences. However, using larger buffers and/or higher resoluted raster images, these discrepancies become larger and larger.



Answer



By making a buffer you approximate the distances (as you do not get a true circle, but a polygonal approximation of one). You can improve the rgeos answer by increasing the number of segments of the circle with the "quadsegs" argument:



nbg.bff2 <- gBuffer(nuremberg, quadsegs=50, width = 20000)
table( unlist(extract(germany, nbg.bff2)) )

0 1
118 91

No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...