Search code examples
rrasterscatter-plottiffr-raster

Computing R2 values according to land cover types between 2 rasters in R?

I have three rasters in r

> lpjre
class      : RasterLayer 
dimensions : 2803, 5303, 14864309  (nrow, ncol, ncell)
resolution : 0.008333333, 0.008333333  (x, y)
extent     : 60.85, 105.0417, 15.95833, 39.31667  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
source     : memory
names      : xxx 
values     : 0, 21  (min, max)

> gcre
class      : RasterLayer 
dimensions : 2803, 5303, 14864309  (nrow, ncol, ncell)
resolution : 0.008333333, 0.008333333  (x, y)
extent     : 60.85, 105.0417, 15.95833, 39.31667  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
source     : memory
names      : layer 
values     : 0, 39.72  (min, max)

and a landcover raster

> tif4
class      : RasterLayer 
dimensions : 2803, 5303, 14864309  (nrow, ncol, ncell)
resolution : 0.008333333, 0.008333333  (x, y)
extent     : 60.85, 105.0417, 15.95833, 39.31667  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
source     : C:/Users/XXXX/landusemaskedme.tif 
names      : landusemaskedme 
values     : 1, 12  (min, max)
attributes :
       ID                          zn
 from:  1 evergreen needleleaf forest
  to : 12                   croplands

I plot a scatterplot between lpjre and gcre according to different land cover classes by

plot(lpjre[tif4==2],gcre[tif4==2])

I would like to find how can I compute the r2 value between lpjre and gcre according to land cover types?

I try this code and gives an error:

> cor(values(lpjre)[tif4==1], values(gcre)[tif4==1], use="complete.obs", method = 'pearson')
Error in values(gcre)[tif4 == 1] : invalid subscript type 'S4'

Reproducible rasters:

library(raster)
ras1 <- raster(matrix(c(1,1,1,2,2,2)))
ras2 <- raster(matrix(c(1,1,1,2,2,2)))

#Generating landcover example data
raster2 <- raster(matrix(c(1,1,1,2,2,2,3,3,3),ncol =3))
raster2 <- as.factor(raster2)

rat <- levels(raster2 )[[1]]
rat[["landcover"]] <- c("land","ocean/lake", "rivers")
levels(raster2 ) <- rat
Solution

  • You can use vector conversion for that. First a reproducible example:

    library(raster)
set.seed(42)

ras1 <- raster(nrow = 10, ncol = 10)
ras2 <- raster(nrow = 10, ncol = 10)
lcc <- raster(nrow = 10, ncol = 10)
ras1[] <- runif(100)
ras2[] <- runif(100)
lcc[] <- sample(c(1,2), 100, replace = TRUE)

    Now you can convert the raster values into vectors by using []. From there on you can use logical selection even between the rasters assuming they have corresponding geometric features. 

    cor(ras1[lcc[] == 2], ras2[lcc[] == 2], use = "complete.obs", method = "pearson")
# [1] -0.1644459