从RasterBrick的多个图层中提取r中每个GPS位置的值

Question

我试图找出USDA CropScape在该位置的日历年、前一个日历年和下一个日历年的作物价值(例如，如果GPS位置的日期是10/1/2016，我想提取该位置2015年、2016年和2017年的作物价值)。

下面是我用来准备GPS和CropScape数据的代码:

library(rgdal)
library(raster)
library(sf)
library(sp)
library(rgeos)
library(adehabitatHR)
library(lubridate)
head(data)
proj <- 32612
data <- st_as_sf(data, coords=c("easting","northing"), dim="XY", crs=proj)
cropScape_2014 <- raster("./data/_raw/cropScape/cropScape_2014.tif")
cropScape_2015 <- raster("./data/_raw/cropScape/cropScape_2015.tif")
cropScape_2016 <- raster("./data/_raw/cropScape/cropScape_2016.tif")
cropScape_2017 <- raster("./data/_raw/cropScape/cropScape_2017.tif")
cropScape_2018 <- raster("./data/_raw/cropScape/cropScape_2018.tif")
cropScape_2019 <- raster("./data/_raw/cropScape/cropScape_2019.tif")
cropScape_2020 <- raster("./data/_raw/cropScape/cropScape_2020.tif")
cropScape_2021 <- raster("./data/_raw/cropScape/cropScape_2021.tif")
cropScape <- stack(cropScape_2014, cropScape_2015, cropScape_2016,
cropScape_2017, cropScape_2018, cropScape_2019,
cropScape_2020, cropScape_2021)
data <- as(data, "Spatial")
popMCP <- mcp(data, percent = 100)
popMCPbuf <- raster::buffer(popMCP, 1000)
cropScape2 <- mask(cropScape, popMCPbuf)
cropScape_cropped <- crop(cropScape2, popMCPbuf)
tz(data$date)
data$calYr <- as.numeric(strftime(data$date, format = "%Y", tz = "MST7MDT"))
dput(data[sample(nrow(data), size = 10),])
dput(head(cropScape_cropped))

这是我在这一点上获取的GPS数据的一小部分:

new("SpatialPointsDataFrame", data = structure(list(id = c("GT74", 
"GT48", "GT52", "GT82", "GT74", "GT52", "GT82", "GT9", "GT45", 
"GT43"), date = structure(c(1607320849, 1544097667, 1588500050, 
1554854505, 1618632063, 1549476050, 1577080837, 1449784850, 1542376850, 
1548309669), class = c("POSIXct", "POSIXt"), tzone = "MST7MDT"), 
id_bioYr = c("GT74_2020", "GT48_2018", "GT52_2019", "GT82_2018", 
"GT74_2020", "GT52_2018", "GT82_2019", "GT9_2015", "GT45_2018", 
"GT43_2018"), state = c(1L, 1L, 3L, 2L, 1L, 1L, 2L, 2L, 2L, 
1L), used = c(0, 0, 0, 1, 0, 0, 0, 0, 0, 0), landcover = c("agricultural", 
"agricultural", "agricultural", "agricultural", "agricultural", 
"shadow", "water", "xericGrass", "agricultural", "agricultural"
), calYr = c(2020, 2018, 2020, 2019, 2021, 2019, 2019, 2015, 
2018, 2019)), row.names = c(471519L, 253950L, 310790L, 49469L, 
479572L, 293620L, 573180L, 627720L, 238551L, 151426L), class = "data.frame"), 
coords.nrs = numeric(0), coords = structure(c(483309.821132995, 
480304.323982496, 479889.104783906, 469493, 458709.035703965, 
475920.645473197, 470714.848866241, 474984.326637862, 473390.631355455, 
479941.915819768, 4866298.15517148, 4867412.65105531, 4865629.72751859, 
4865043, 4864209.38410713, 4863155.30805073, 4865581.17021548, 
4865450.06432659, 4865829.2649267, 4861889.90054007), .Dim = c(10L, 
2L), .Dimnames = list(NULL, c("coords.x1", "coords.x2"))), 
bbox = structure(c(458709.035703965, 4861889.90054007, 483309.821132995, 
4867412.65105531), .Dim = c(2L, 2L), .Dimnames = list(c("coords.x1", 
"coords.x2"), c("min", "max"))), proj4string = new("CRS", 
projargs = "+proj=utm +zone=12 +datum=WGS84 +units=m +no_defs"))

和作物数据(不确定这是否有用/人们需要什么?):

structure(c(NA_integer_, NA_integer_, NA_integer_, NA_integer_, 
NA_integer_, NA_integer_, NA_integer_, NA_integer_, NA_integer_, 
NA_integer_, NA_integer_, NA_integer_, NA_integer_, NA_integer_, 
NA_integer_, NA_integer_, NA_integer_, NA_integer_, NA_integer_, 
NA_integer_, NA_integer_, NA_integer_, NA_integer_, NA_integer_, 
NA_integer_, NA_integer_, NA_integer_, NA_integer_, NA_integer_, 
NA_integer_, NA_integer_, NA_integer_, NA_integer_, NA_integer_, 
NA_integer_, NA_integer_, NA_integer_, NA_integer_, NA_integer_, 
NA_integer_, NA_integer_, NA_integer_, NA_integer_, NA_integer_, 
NA_integer_, NA_integer_, NA_integer_, NA_integer_, NA_integer_, 
NA_integer_, NA_integer_, NA_integer_, NA_integer_, NA_integer_, 
NA_integer_, NA_integer_, NA_integer_, NA_integer_, NA_integer_, 
NA_integer_, NA_integer_, NA_integer_, NA_integer_, NA_integer_, 
NA_integer_, NA_integer_, NA_integer_, NA_integer_, NA_integer_, 
NA_integer_, NA_integer_, NA_integer_, NA_integer_, NA_integer_, 
NA_integer_, NA_integer_, NA_integer_, NA_integer_, NA_integer_, 
NA_integer_, NA_integer_, NA_integer_, NA_integer_, NA_integer_, 
NA_integer_, NA_integer_, NA_integer_, NA_integer_, NA_integer_, 
NA_integer_, NA_integer_, NA_integer_, NA_integer_, NA_integer_, 
NA_integer_, NA_integer_, NA_integer_, NA_integer_, NA_integer_, 
NA_integer_, NA_integer_, NA_integer_, NA_integer_, NA_integer_, 
NA_integer_, NA_integer_, NA_integer_, NA_integer_, NA_integer_, 
NA_integer_, NA_integer_, NA_integer_, NA_integer_, NA_integer_, 
NA_integer_, NA_integer_, NA_integer_, NA_integer_, NA_integer_, 
NA_integer_, NA_integer_, NA_integer_, NA_integer_, NA_integer_, 
NA_integer_, NA_integer_, NA_integer_, NA_integer_, NA_integer_, 
NA_integer_, NA_integer_, NA_integer_, NA_integer_, NA_integer_, 
NA_integer_, NA_integer_, NA_integer_, NA_integer_, NA_integer_, 
NA_integer_, NA_integer_, NA_integer_, NA_integer_, NA_integer_, 
NA_integer_, NA_integer_, NA_integer_, NA_integer_, NA_integer_, 
NA_integer_, NA_integer_, NA_integer_, NA_integer_, NA_integer_, 
NA_integer_, NA_integer_, NA_integer_, NA_integer_, NA_integer_, 
NA_integer_), .Dim = c(20L, 8L), .Dimnames = list(NULL, c("cropScape_2014", 
"cropScape_2015", "cropScape_2016", "cropScape_2017", "cropScape_2018", 
"cropScape_2019", "cropScape_2020", "cropScape_2021")))

我一直在尝试编写一个循环，将我的三个值提取到数据中的新列中，但我真的不知道如何开始。

理想情况下，在循环结束时，将有3个新列添加到数据中("crop_calYr-1"， " crop_caly"， "crop_calYr+1")，并且它们将从" cropscape_裁剪。"的正确层填充。我只知道我得加上" NA";crop_calYr + 1"当GPS位置从2021年开始

Answer 1

我无法重新创建您的栅格数据。我刚下载了弗里蒙特县几年来的数据。因为看起来你的点就在那里。我在这里使用的是raster包，但我更愿意使用terra。

library(raster)
l1 <- list.files("Downloads/polygonclip", pattern="tif$", full.names
=T)
cd1 <- raster::stack(l1)
# Should be able to change the categories in the stack
# so that the extract returns character value, e.g. 
# "barley" instead of 152 with something like this. 
# Get levels
ca1 <- levels(cd1) 
# MOdify to get desired factors
ca2 <- ca1[[1]][,c(1,5)]
# Assign
levels(cd2, layer=c(1:4)) <- ca2
# ... But that way doesn't work for me for some reason, 
# so I did it this way. You'll need to do this for each
# year/layer in the stack
levels(cd1[[1]]) <- ca2
levels(cd1[[2]]) <- ca2
levels(cd1[[3]]) <- ca2
levels(cd1[[4]]) <- ca2
# Recreate the points data shared
pts1 <- new("SpatialPointsDataFrame", data = structure(list(id = c("GT74", 
"GT48", "GT52", "GT82", "GT74", "GT52", "GT82", "GT9", "GT45", 
"GT43"), date = structure(c(1607320849, 1544097667, 1588500050, 
1554854505, 1618632063, 1549476050, 1577080837, 1449784850, 1542376850, 
1548309669), class = c("POSIXct", "POSIXt"), tzone = "MST7MDT"), 
id_bioYr = c("GT74_2020", "GT48_2018", "GT52_2019", "GT82_2018", 
"GT74_2020", "GT52_2018", "GT82_2019", "GT9_2015", "GT45_2018", 
"GT43_2018"), state = c(1L, 1L, 3L, 2L, 1L, 1L, 2L, 2L, 2L, 
1L), used = c(0, 0, 0, 1, 0, 0, 0, 0, 0, 0), landcover = c("agricultural", 
"agricultural", "agricultural", "agricultural", "agricultural", 
"shadow", "water", "xericGrass", "agricultural", "agricultural"
), calYr = c(2020, 2018, 2020, 2019, 2021, 2019, 2019, 2015, 
2018, 2019)), row.names = c(471519L, 253950L, 310790L, 49469L, 
479572L, 293620L, 573180L, 627720L, 238551L, 151426L), class = "data.frame"), 
coords.nrs = numeric(0), coords = structure(c(483309.821132995, 
480304.323982496, 479889.104783906, 469493, 458709.035703965, 
475920.645473197, 470714.848866241, 474984.326637862, 473390.631355455, 
479941.915819768, 4866298.15517148, 4867412.65105531, 4865629.72751859, 
4865043, 4864209.38410713, 4863155.30805073, 4865581.17021548, 
4865450.06432659, 4865829.2649267, 4861889.90054007), .Dim = c(10L, 
2L), .Dimnames = list(NULL, c("coords.x1", "coords.x2"))), 
bbox = structure(c(458709.035703965, 4861889.90054007, 483309.821132995, 
4867412.65105531), .Dim = c(2L, 2L), .Dimnames = list(c("coords.x1", 
"coords.x2"), c("min", "max"))), proj4string = new("CRS", 
projargs = "+proj=utm +zone=12 +datum=WGS84 +units=m +no_defs"))
# Removed one of the points because it was for the year 2015 
# which threw everything off. 
pts1 <- pts1[-8,]
# Need same projection for extract        
pts1 <- spTransform(pts1, crs(cd1))
# Set Z dimension (years) for matching in next function
# Note the years used
cd1 <- setZ(cd1, 2018:2021, "years")
# This takes a year and returns the indices of that year
# from the list of years. The indices are for the year
# prior to the selected year, the year selected, and the
# year after selected year.
# Note the years used. 2018 was my oldest, but I think
# your data goes to 2015, so you'll have to change each
# 2018
rysel <- function(yind) {
if(yind==2018) {
yex = c(1, 2)
} else if(yind==2021) {
yex = c(3, 4)
} else if(2018 < yind & yind < 2021) {
yex = c(yind-2018, yind-2018+1, yind-2018+2)
}
return(yex)
}
# This takes the index of the row in the points data
# and uses the rysel function to subset the stack and 
# then extract.
extryr <- function(x) {
pty = pts1$calYr[x]
ytmp = rysel(pty)
yextr = extract(cd1[[ytmp]], pts1[x,])
# Pad the extract for either end
# Note the years used
if(pty==2018){
yextr = c(yextr, "NA")
} else if(pty==2021) {
yextr = c("NA",yextr)
}
return(yextr)
}
# Data frame with column labels for holding extract values
df2 <- data.frame(calYrprior = NA, calYrcur = NA, calYrnext = NA)
# 
for(i in 1:nrow(pts1)) df2[i,] = extryr(i)
# Convert the values in df2 to numeric
# You can ignore the warning
df2 <- sapply(df2, as.numeric)
Warning messages:
1: In extryr(i) : NAs introduced by coercion
2: In extryr(i) : NAs introduced by coercion
3: In extryr(i) : NAs introduced by coercion
> df2
calYrprior calYrcur calYrnext
[1,]        152      152       152
[2,]         21       21        NA
[3,]        152      152       152
[4,]        152      152        21
[5,]         NA       23        43
[6,]          0        0         0
[7,]          0        0         0
[8,]        152      152        NA
[9,]         NA       NA        NA

# Combine the pts1 data and the extract data
pts1 <- cbind(pts1, df2)

你需要把我标注的年份改一下。如果在栅格和点上的年份有差异，这将不起作用。例如，如果堆栈有2014年，2016年至2019年，而点数涵盖2014年至2019年，则必须修改此内容。

提取仍然返回数值而不是字符/因子值。但是您可以根据需要执行替换来更改这些。

Answer 2

正如John Polo所建议的那样，使用"terra"要容易得多。(替换"栅格")

示例数据

library(terra)
f <- system.file("ex/elev.tif", package="terra")
r <- rast(f)
x <- c(r*1, r*2, r*3, r*4, r*5)
time(x, "years") <- 2001:2005
names(x) <- paste0("Y", 2001:2005)
d <- data.frame(lon=c(5.867235,  6.319561), 
lat=c(49.85114, 49.55199), year=c(2003, 2004))

解决方案

# match the year in your data.frame with that of the SpatRaster layers
i <- match(d$year, time(x))
## or 
# i <- match(d$year, 2001:2005)

如果您只需要匹配年份的值，您可以执行

ey <- extract(x, d[, c("lon","lat")], layer=i)
ey
#  ID layer value
#1  1 Y2003  1473
#2  2 Y2004   692

还可以得到前后年份:

## extract values for all years
e <- extract(x, d[, c("lon","lat")], ID=FALSE)
## find the layers (years) for each point
j <- cbind(rep(1:length(d$year), each=3), 
rep(i, each=3) + c(-1:1))
## subset the extracted data
e[j]
#[1]  982 1473 1964  519  692  865
## to keep track of what is what
b <- data.frame(point=j[,1], year=time(x)[j[,2]], value=e[j])
b
#  point year value
#1     1 2002   982
#2     1 2003  1473
#3     1 2004  1964
#4     2 2003   519
#5     2 2004   692
#6     2 2005   865
reshape(b, direction="wide", idvar="point", timevar="year")
#  point value.2002 value.2003 value.2004 value.2005
#1     1        982       1473       1964         NA
#4     2         NA        519        692        865 

这是几天前的一个类似(但更复杂)的问题。