如何为每年创建一层的 RasterStack答案

【问题标题】：How to create RasterStack with one layer for each year如何为每年创建一层的 RasterStack
【发布时间】：2019-07-17 08:39:32
【问题描述】：

我有一个 sf 对象，其中包含对多个单位的变量的年度观察。我想对此进行栅格化，以便从 1946 年到 2005 年每年都有一个 RasterLayer （观察的开始年份不同）。数据如下所示：

> head(data)
Simple feature collection with 6 features and 4 fields
geometry type:  POINT
dimension:      XY
bbox:           xmin: 0.5 ymin: 26.63 xmax: 0.5 ymax: 26.63
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs
# A tibble: 6 x 5
   id   year       var1        geometry
  <dbl> <int>     <dbl>       <POINT [°]>
1   615  1957         1        (0.5 26.63)
2   615  1958         1        (0.5 26.63)
3   615  1959         1        (0.5 26.63)
4   615  1960         1        (0.5 26.63)
5   615  1961         1        (0.5 26.63)
6   615  1962         1        (0.5 26.63)

每个RasterLayer 都应包含所有几何点的var1 是否存在的信息。

我尝试使用以下代码，但这会创建一个两层 RasterBrick，而不是我想要的 59 层 RasterStack。

r <- raster::rasterize(data, raster, field = c("var1", "year"), fun = "first")

【问题讨论】：

您可以使用 for 循环对每一年进行栅格化，然后在之后堆叠栅格图层。
该功能未设置为您可以循环，或将您的点更改为宽格式，其中易年是一个变量。

标签： r r-raster sf rasterizing

【解决方案1】：

我使用循环解决了它，需要一些时间来运行 40 000 个位置，但它可以工作

setwd("E:/MA_2Try/annual_NDVI/average")
grids <- list.files("E:/MA_2Try/annual_NDVI/average" , pattern = "*.tif$",full.names=TRUE)
NDVI_stack=stack(grids)

#Add dates to stack layers
dt<-as.data.frame(
  as.Date(c('2002-01-01','2003-01-01','2004-01-01','2005-01-01','2006-01-01',
            '2007-01-01','2008-01-01','2009-01-01','2010-01-01','2011-01-01',
            '2012-01-01','2013-01-01','2014-01-01','2015-01-01','2016-01-01',
            '2017-01-01','2018-01-01','2019-01-01')))
years <- format(as.Date(dt[,1], format="%d/%m/%Y"),"%Y")
NDVI_stack <- setZ(NDVI_stack, years,"year")

#extract NDVI values
ex <- list()

for(i in 1:nrow(data)) {
  
  if (data$year[i] == 2009) {
    NDVI <- raster::extract(NDVI_stack[[8]], data[, 5:6])
  } 
  else if (data$year[i]==2010) {
    NDVI <- raster::extract(NDVI_stack[[9]], data[, 5:6])
  } 
  else if (data$year[i]==2011) {
    NDVI <- raster::extract(NDVI_stack[[10]], data[, 5:6])
  } 
  else if (data$year[i]==2012) {
    NDVI <- raster::extract(NDVI_stack[[11]], data[, 5:6])
  } 
  else if (data$year[i]==2013) {
    NDVI <- raster::extract(NDVI_stack[[12]], data[, 5:6])
  } 
  else if (data$year[i]==2014) {
    NDVI <- raster::extract(NDVI_stack[[13]], data[, 5:6])
  }  
  else if (data$year[i]==2015) {
    NDVI <- raster::extract(NDVI_stack[[14]], data[, 5:6])
  }
  else if (data$year[i]==2016) {
    NDVI <- raster::extract(NDVI_stack[[15]], data[, 5:6])
  }
  else if (data$year[i]==2017) {
    NDVI <- raster::extract(NDVI_stack[[16]], data[, 5:6])
  }
  else if (data$year[i]==2018) {
    NDVI <- raster::extract(NDVI_stack[[17]], data[, 5:6])
  }
  else if (data$year[i]==2019) {
    NDVI <- raster::extract(NDVI_stack[[18]], data[, 5:6])
  }
  else if (data$year[i]==2020) {
    NDVI <- raster::extract(NDVI_stack[[18]], data[, 5:6])
  }
  else if (data$year[i]==2021) {
    NDVI <- raster::extract(NDVI_stack[[18]], data[, 5:6])
  }
  ex[[i]] <- list(id=data$id[i], y=data$y_[i], NDVI= NDVI[i])
  
  progress(i, 41290)
  Sys.sleep(0.02)
  if (i == 41290) message("Done!")
  
}
output <- do.call(rbind.data.frame, ex)
setwd("E:/Masterthesis/wildcat_data/second_order/seasons")
write.csv(output,'extracted.NDVI.summer.csv')

【讨论】：