如何从R中的栅格堆栈中提取时间序列答案

【问题标题】：How to extract time series from a raster stack in R如何从R中的栅格堆栈中提取时间序列
【发布时间】：2019-07-06 22:15:58
【问题描述】：

我已经从 MODIS 数据的 NDVI 层创建了一个 RasterStack。现在我想从这些数据的不同位置提取时间序列数据，以便我可以使用 BFAST/greenbrown 包来估计趋势和断点。这是我创建堆栈的方式：

#runGdal(Job="testJob","MYD13Q1",begin = "2018.01.09", end = "2018.12.27",
#        tileH = 26:29, tileV = 4:7
#        , SDSstring = "1000000000000000000000") 

###NDVI files path
NDVI_files_path <- "/media/MyData/Data/MODIS/PROCESSED/MYD13Q1.006_20190527193158"
all_NDVI_files <- list.files(NDVI_files_path,
                            full.names = TRUE,
                            pattern = ".tif$")
all_NDVI_files

### Raster Stack
NDVI_stack <- stack(all_NDVI_files)

如何提取栅格堆栈中任何特定区域的时间序列数据？

【问题讨论】：

minimal reproducible example 将包含任何需要的 library 调用和一个可用于测试的小示例。
我会看看函数 raster::extract()

标签： r time-series r-raster

【解决方案1】：

您可以使用lubridate 和rts 创建一个RasterStackTS 对象，如下所示：

#Load libraries
library(rts)
library(lubridate)

#Reproducible example (use your files here)
all_NDVI_files = c('MOD13Q1.006__250m_16_days_EVI_doy2000049_aid0001.tif',
               'MOD13Q1.006__250m_16_days_EVI_doy2000065_aid0001.tif',
               'MOD13Q1.006__250m_16_days_EVI_doy2000081_aid0001.tif',
               'MOD13Q1.006__250m_16_days_EVI_doy2000097_aid0001.tif',
               'MOD13Q1.006__250m_16_days_EVI_doy2000113_aid0001.tif',
               'MOD13Q1.006__250m_16_days_EVI_doy2000129_aid0001.tif',
               'MOD13Q1.006__250m_16_days_EVI_doy2000145_aid0001.tif',
               'MOD13Q1.006__250m_16_days_EVI_doy2000161_aid0001.tif',
               'MOD13Q1.006__250m_16_days_EVI_doy2000177_aid0001.tif',
               'MOD13Q1.006__250m_16_days_EVI_doy2000193_aid0001.tif',
               'MOD13Q1.006__250m_16_days_EVI_doy2000209_aid0001.tif',
               'MOD13Q1.006__250m_16_days_EVI_doy2000225_aid0001.tif',
               'MOD13Q1.006__250m_16_days_EVI_doy2000241_aid0001.tif',
               'MOD13Q1.006__250m_16_days_EVI_doy2000257_aid0001.tif',
               'MOD13Q1.006__250m_16_days_EVI_doy2000273_aid0001.tif',
               'MOD13Q1.006__250m_16_days_EVI_doy2000289_aid0001.tif',
               'MOD13Q1.006__250m_16_days_EVI_doy2000305_aid0001.tif',
               'MOD13Q1.006__250m_16_days_EVI_doy2000321_aid0001.tif',
               'MOD13Q1.006__250m_16_days_EVI_doy2000337_aid0001.tif',
               'MOD13Q1.006__250m_16_days_EVI_doy2000353_aid0001.tif',
               'MOD13Q1.006__250m_16_days_EVI_doy2001001_aid0001.tif',
               'MOD13Q1.006__250m_16_days_EVI_doy2001017_aid0001.tif',
               'MOD13Q1.006__250m_16_days_EVI_doy2001033_aid0001.tif',
               'MOD13Q1.006__250m_16_days_EVI_doy2001049_aid0001.tif',
               'MOD13Q1.006__250m_16_days_EVI_doy2001065_aid0001.tif',
               'MOD13Q1.006__250m_16_days_EVI_doy2001081_aid0001.tif',
               'MOD13Q1.006__250m_16_days_EVI_doy2001097_aid0001.tif',
               'MOD13Q1.006__250m_16_days_EVI_doy2001113_aid0001.tif',
               'MOD13Q1.006__250m_16_days_EVI_doy2001129_aid0001.tif',
               'MOD13Q1.006__250m_16_days_EVI_doy2001145_aid0001.tif',
               'MOD13Q1.006__250m_16_days_EVI_doy2001161_aid0001.tif',
               'MOD13Q1.006__250m_16_days_EVI_doy2001177_aid0001.tif',
               'MOD13Q1.006__250m_16_days_EVI_doy2001193_aid0001.tif',
               'MOD13Q1.006__250m_16_days_EVI_doy2001209_aid0001.tif',
               'MOD13Q1.006__250m_16_days_EVI_doy2001225_aid0001.tif',
               'MOD13Q1.006__250m_16_days_EVI_doy2001241_aid0001.tif',
               'MOD13Q1.006__250m_16_days_EVI_doy2001257_aid0001.tif',
               'MOD13Q1.006__250m_16_days_EVI_doy2001273_aid0001.tif',
               'MOD13Q1.006__250m_16_days_EVI_doy2001289_aid0001.tif',
               'MOD13Q1.006__250m_16_days_EVI_doy2001305_aid0001.tif')

#Depending on file format, extract dates (this example uses MODIS 16-day composite NDVI)
ndvi.time = data.frame(year=substr(basename(all_NDVI_files),34,37),
                julD=substr(basename(all_NDVI_files),38,40))
ndvi.time$dateJ = paste(ndvi.time$year,ndvi.time$julD,sep='-')
ndvi.time$julD = parse_date_time(ndvi.time$dateJ,'y-j')

# create your RasterStackTS object
# Use your actual rasterstack here i.e. it's not reproducible with above code
ndvi.rts = rts(NDVI_stack ,ndvi.time$julD) 
ndvi.rts
plot(ndvi.rts)

希望这会有所帮助！

【讨论】：

感谢您的示例。我创建了时间序列对象，但是当我想在bfastmonitor 或greenbrown 中使用它时。它给了我ndvi.rts should be from class 'ts' 的错误。如何将创建的rts 对象转换为ts 对象，以便绘制两个包文档中显示的趋势图？
据我所知greenbrown 不使用RasterStackTS 对象。只需使用原始的RasterStack 或RasterBrick 并添加开始和频率变量，如下所示：TrendRaster(NDVI_stack, start=min(ndvi.time$julD), freq=23, method='AAT', breaks=4)。 16 天复合材料的频率为 23（每年）。