If buffer is in meters the code works, but if I use degrees, it won't.
Indeed; the buffer cannot be specified in degrees. ?extract, buffer says: "If the data are not projected (latitude/longitude), the unit should be meters. Otherwise it should be in map-units (typically also meters)."
Do you really want the buffer to be in degrees?
Is so, temporarily set the CRS to UTM, do the extract, and set it back.
cr <- crs(r)
crs(r) <- '+proj=utm +zone=10 +datum=WGS84'
e <- extract(r, xy, buffer=0.5)
crs(r) <- cr
I tried to "parallelize" the function rasterize using the R package parallel in this way:
- split the SpatialPolygonsDataFrame object in
n parts
rasterize every part separately- merge all the parts into one raster
In my computer, the parallelized rasterize function took 2.75 times less than the no-parallelized rasterize function.
Note: the code below download a polygon shapefile (~26.2 MB) from the web. You can use any SpatialPolygonDataFrame object. This is only an example.
Load libraries and example data:
# Load libraries
library('raster')
library('rgdal')
# Load a SpatialPolygonsDataFrame example
# Load Brazil administrative level 2 shapefile
BRA_adm2 <- raster::getData(country = "BRA", level = 2)
# Convert NAMES level 2 to factor
BRA_adm2$NAME_2 <- as.factor(BRA_adm2$NAME_2)
# Plot BRA_adm2
plot(BRA_adm2)
box()
# Define RasterLayer object
r.raster <- raster()
# Define raster extent
extent(r.raster) <- extent(BRA_adm2)
# Define pixel size
res(r.raster) <- 0.1
Figure 1: Brazil SpatialPolygonsDataFrame plot
Simple thread example
# Simple thread -----------------------------------------------------------
# Rasterize
system.time(BRA_adm2.r <- rasterize(BRA_adm2, r.raster, 'NAME_2'))
Time in my laptop:
# Output:
# user system elapsed
# 23.883 0.010 23.891
Multithread thread example
# Multithread -------------------------------------------------------------
# Load 'parallel' package for support Parallel computation in R
library('parallel')
# Calculate the number of cores
no_cores <- detectCores() - 1
# Number of polygons features in SPDF
features <- 1:nrow(BRA_adm2[,])
# Split features in n parts
n <- 50
parts <- split(features, cut(features, n))
# Initiate cluster (after loading all the necessary object to R environment: BRA_adm2, parts, r.raster, n)
cl <- makeCluster(no_cores, type = "FORK")
print(cl)
# Parallelize rasterize function
system.time(rParts <- parLapply(cl = cl, X = 1:n, fun = function(x) rasterize(BRA_adm2[parts[[x]],], r.raster, 'NAME_2')))
# Finish
stopCluster(cl)
# Merge all raster parts
rMerge <- do.call(merge, rParts)
# Plot raster
plot(rMerge)
Figure 2: Brazil Raster plot
Time in my laptop:
# Output:
# user system elapsed
# 0.203 0.033 8.688
More info about parallelization in R:
Best Answer
Your data is in Lat-long (EPSG:4326 by the look of it) so you are making a buffer of 1000 degrees which is going to wrap round the earth about three times!
You need to transform your spatial object to a different coordinate system based on metres (or km, or whatever units your 1000 is) using
spTransform. I suspect for France there's a reasonable UTM zone with corresponding EPSG code. That will give you metres for a 1km zone with width=1000.Note that for a whole country like France the outline may be very detailed and you might want to use gSimplify first to get a simpler outline.