Raster in R – Synchronizing RasterStack and SpatialPolygonsDataFrame to the Same Coordinate System

rrasterrgdal

I have a shapefile of 8094 Italian municipalities which I imported into R as a SpatialPolygonsDataFrame using readOGR. When I plot this, I get the following map:

I also have a RasterStack of global monthly temperatures.

I want to overlay the rasterstack with the shapefile to extract the values for each municipality. I have done this many times before (with other shapefiles) using the extract command in R, i.e., monthly_temperatures <- raster::extract(tmp_month, Italy_map)

However, this time it doesn't work because the RasterStack and the SpatialPolygonsDataFrame use different coordinate systems. The proj4string of the Italy map is +proj=utm +zone=32 +datum=WGS84 +units=m +no_defs and the bbox is as follows:

    min     max
x  313279.3 1312016
y 3933520.1 5220291

Of the RasterStack the proj4string is +proj=longlat +datum=WGS84 +no_defs and the bbox is as follows:

    min     max
x  -180    180
y   -90     90

They clearly do not overlap. Does someone know how to fix this issue? I have looked for different shapefiles of Italian municipalities online but the one from ISTAT is the only one I have found. I am working in RStudio.

Best Answer

As far I know, R doesn't project on the fly while processing with multiple layers. All layers should share CRS for any process.

For projecting CRS, use spTransform. Always project the vector file over the raster file:

Italy_mapWGS84 <- spTransform(Italy_map, CRS('+proj=longlat +datum=WGS84 +no_defs'))

monthly_temperatures <- raster::extract(tmp_month, Italy_mapWGS84)

Related Solutions

[GIS] Extract values from raster using point data with buffer in R

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

Process Vector to Raster Faster with R – Speed Optimization Techniques

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