MODISTools的時間序列

Question

我需要獲取完整的EVI時間序列以及日期和質量信息。在執行MODISSubsets()之後，原始數據可用，但沒有像MODISSummaries()那樣以相當好的方式處理。

MODISSummaries()但是，考慮到質量信息，縮短了彙總統計的時間序列。

有沒有一種方法可以從原始數據中提取每個圖塊的時間序列（請參閱下面的數據框crude）？如果這可能會返回一個數據幀列表，其中每個數據幀代表一個塊並保存EVI（或其他變量）的數據，其日期和質量標誌，那將是非常好的。

具體來說，後做以下...

savedir <- './' 

modis.subset <- data.frame( 
    lat  = 11.3175, 
    long  = 47.1167, 
    end.date = "2016-09-29" 
) 

MODISSubsets(
    LoadDat = modis.subset, 
    Products = "MOD13Q1", 
    Bands  = c("250m_16_days_EVI", "250m_16_days_pixel_reliability"), 
    Size  = c(1,1), 
    StartDate = FALSE, 
    SaveDir = savedir, 
    TimeSeriesLength = 3 
) 

crude <- read.csv("./Lat47.11670Lon11.31750Start2013-01-01End2016-09-29___MOD13Q1.asc", header = FALSE, as.is = TRUE) 

...你將如何去像

nice <- list(lonX1_latY1=data.frame(date=..., var=..., qual=...), lonX2_latX2=... ) 

...？

Answer 1

總之，我錯過了ExtractTile()將可用與返回值MODISTimeSeries()。我的解決方法是基於使用ExtractTile()結合讀取ASCII文件的輸出。這裏是我爲我的目的而工作，返回一個列表，其中包含包含所有下載的MODIS數據的矩陣（npixels_lon，npixels_lat，n_timesteps），在這種情況下爲EVI;包含像素可靠性代碼的相同維度的矩陣;和長度的矢量n_timesteps保持中心像素信息，如果其質量標誌是0或其周圍像素的平均值，否則：

read_crude_modis <- function(filn, savedir, expand_x, expand_y){ 

    # arguments: 
    # filn: file name of ASCII file holding MODIS "crude" data 
    # savedir: directory, where to look for that file 
    # expand_x : number of pixels to the right and left of centre 
    # expand_y : number of pixels to the top and bottom of centre 

    # MODIS quality flags: 
    # -1 Fill/No Data Not Processed 
    # 0 Good Data Use with confidence 
    # 1 Marginal data Useful, but look at other QA information 
    # 2 Snow/Ice Target covered with snow/ice 
    # 3 Cloudy Target not visible, covered with cloud 

    library(MODISTools) 

    ScaleFactor <- 0.0001 # applied to output variable 
    ndayyear <- 365 

    ## Read dowloaded ASCII file 
    crude <- read.csv(paste(savedir, filn, sep=""), header = FALSE, as.is = TRUE) 
    crude <- rename(crude, c("V1"="nrows", "V2"="ncols", "V3"="modislon_ll", "V4"="modislat_ll", "V5"="dxy_m", "V6"="id", "V7"="MODISprod", "V8"="yeardoy", "V9"="coord", "V10"="MODISprocessdatetime")) 

    ## this is just read to get length of time series and dates 
    tseries <- MODISTimeSeries(savedir, Band = "250m_16_days_EVI") 
    ntsteps <- dim(tseries[[1]])[1] 
    tmp  <- rownames(tseries[[1]]) 
    time  <- data.frame(yr=as.numeric(substr(tmp, start=2, stop=5)), doy=as.numeric(substr(tmp, start=6, stop=8))) 
    time$dates <- as.POSIXlt(as.Date(paste(as.character(time$yr), "-01-01", sep="")) + time$doy - 1) 
    time$yr_dec<- time$yr + (time$doy - 1)/ndayyear 

    ## get number of products for which data is in ascii file (not used) 
    nprod <- dim(crude)[1]/ntsteps 
    if ((dim(crude)[1]/nprod)!=ntsteps) { print("problem") } 

    ## re-arrange data 
    if (dim(crude)[2]==11 && expand_x==0 && expand_y==0){ 
    ## only one pixel downloaded 
    nice_all  <- as.matrix(crude$V11[1:ntsteps], dim(1,1,ntsteps)) * ScaleFactor  ## EVI data 
    nice_qual_flg <- as.matrix(crude$V11[(ntsteps+1):(2*ntsteps)], dim(1,1,ntsteps))  ## pixel reliability data 

    } else if (dim(crude)[2]>11){ 
    ## multiple pixels downloaded 
    # nice <- ExtractTile(Data = tseries, Rows = c(crude$nrows,expand_y), Cols = c(crude$ncols,expand_x), Grid = TRUE) ## > is not working: applying ExtractTile to return of MODISTimeSeries 
    nice_all  <- ExtractTile(Data = crude[1:ntsteps,11:dim(crude)[2]] * ScaleFactor, Rows = c(crude$nrows[1],expand_y), Cols = c(crude$ncols[1],expand_x), Grid = TRUE) 
    nice_qual_flg <- ExtractTile(Data = crude[(ntsteps+1):(2*ntsteps),11:dim(crude)[2]], Rows = c(crude$nrows[1],expand_y), Cols = c(crude$ncols[1],expand_x), Grid = TRUE) 

    } else { 

    print("Not sufficient data downloaded. Adjust expand_x and expand_y.") 

    } 

    ## Clean data for centre pixel: in case quality flag is not '0', use mean of all 8 surrounding pixels 
    if (expand_x==1 && expand_y==1){ 
    nice_centre <- nice_all[2,2,] 
    nice_centre[ which(nice_qual_flg[2,2,]!=0) ] <- apply(nice_all[,,which(nice_qual_flg[2,2,]!=0)], c(3), FUN=mean) 
    } 

    modis <- list(nice_all=nice_all, nice_centre=nice_centre, nice_qual_flg=nice_qual_flg, time=time) 
    return(modis) 

}