Как извлечь временные ряды из растрового стека в R

Question

Я создал RasterStack из слоев NDVI данных MODIS. Теперь я хочу извлечь данные временных рядов из разных мест этих данных, чтобы я мог использовать пакет 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)

Как извлечь данные временных рядов для любой конкретной области в стеке растров?

Answer 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)

Надеюсь, это поможет!