Extract function in R giving different results for stacked and individual raster

user_9812021 0
Viewed 152

I am getting some differences in extract function in R. Could anyone tell me how to rectify it? In the example below, I have used stacked raster and individual raster for extraction. But I get NA values for some points in individual raster however, stacked raster has no NA values.

The files are in this link: https://drive.google.com/drive/folders/1lmR5YVMw56vTV51855EAWB9htjymZRsT?usp=sharing

setwd("C:/Users/Dell/Desktop/Extract")

library (raster)
b4clip <- raster("wc2.1_30sec_GRD_bio_4.tif")
s6 <- raster("RiverATLAS_30sec_GRD_hft_ix_c09.tif")
s<-stack(b4clip,s6)

library(readxl)
library(sf)
DataSpecies = read_excel("new_extract.xlsx", sheet = "all",
                         col_types = c("numeric"))  
DataSpecies <- DataSpecies [,c("X","Y")]

#Extract for stack
Stacked<-extract(s, DataSpecies) #uncheck tidyr and tidytext while doing this
(Stacked_xy<-cbind(Stacked,DataSpecies))
table(is.na(Stacked_xy))

S6ed = extract(s6, DataSpecies)
(S6ed_xy<-cbind(S6ed,DataSpecies))
table(is.na(S6ed_xy)) #shows 2 NAs
r-raster raster extract r
1 Answers

The apparent inconsistency in the output values has to do with the order the values are extracted, when the point falls just in between the four cells, as in your case. Apparently, extract takes the lowest value from the neighboring cells; and starts reading them in a certain order. The chosen cell in the first layer will be the one for the second (this may not be a deterministic algorithm). Function manual only says:

If y represents points, extract returns the values of a Raster* object for the cells in which a set of points fall

The following code inverts the order of your stack, yielding NAs. A workaround would be to use a buffer in the extract function, as proposed in the following code.

enter image description here

library(raster)
library(readxl)
library(sf)

b4clip <- raster("wc2.1_30sec_GRD_bio_4.tif")
s6 <- raster("RiverATLAS_30sec_GRD_hft_ix_c09.tif")
s_inverse <- stack(s6, b4clip)
s = stack( b4clip, s6)

DataSpecies = read_excel("new_extract.xlsx", sheet = "all",
                         col_types = c("numeric"))  
coordinates(DataSpecies) = c("X", "Y")
crs(DataSpecies) = crs(s6)

# crop to make more manageable layers
s6 = crop(s6, extent(DataSpecies[c(3,5),]) + 1)
b4clip = crop(b4clip, extent(DataSpecies[c(3,5),]) + 1)

# this outputs no NAs
extract(s, DataSpecies)
     wc2.1_30sec_GRD_bio_4 RiverATLAS_30sec_GRD_hft_ix_c09
[1,]              507.5884                             155
[2,]              452.5105                             178
[3,]              439.9914                             145
[4,]              431.8062                             215
[5,]              428.0585                             100
[6,]              423.0356                             162
[7,]              419.9213                             140

# this outputs NAs
extract(s_inverse, DataSpecies)
     RiverATLAS_30sec_GRD_hft_ix_c09 wc2.1_30sec_GRD_bio_4
[1,]                             155              507.5884
[2,]                             233              452.7849
[3,]                              NA              443.0755
[4,]                             215              433.1516
[5,]                              NA                    NA
[6,]                              93              425.1014
[7,]                             119              422.7209


# No NA
extract(s, DataSpecies, method = "bilinear")
# one NA
extract(s_inverse, DataSpecies, method = "bilinear")

# no more NAs, (buffer units are meters, and the buffered polygon should cross cell center)
extract(s_inverse, DataSpecies, buffer = 1000, fun = mean)
     RiverATLAS_30sec_GRD_hft_ix_c09 wc2.1_30sec_GRD_bio_4
[1,]                        155.0000              506.7434
[2,]                        205.5000              452.6477
[3,]                        146.6667              443.5777
[4,]                        215.0000              432.4789
[5,]                        100.0000              428.1727
[6,]                        139.0000              424.7256
[7,]                        130.0000              422.3499

# apparently it takes the lowest value from the first layer
par(mfrow = c(1,2))
c_lim = coordinates(DataSpecies[3,])
plot(s6,  xlim = c(c_lim[1] - .015, c_lim[1] + .015), ylim = c(c_lim[2] - .01, c_lim[2] + .01))
plot(DataSpecies[3,], add = T)
text(s6,  xlim = c(c_lim[1] - .015, c_lim[1] + .015), ylim = c(c_lim[2] - .01, c_lim[2] + .01))

plot(b4clip,  xlim = c(c_lim[1] - .015, c_lim[1] + .015), ylim = c(c_lim[2] - .01, c_lim[2] + .01))
plot(DataSpecies[3,], add = T)
text(b4clip,  xlim = c(c_lim[1] - .015, c_lim[1] + .015), ylim = c(c_lim[2] - .01, c_lim[2] + .01))
user_82559550
Related