8-дневный MODIS растр до среднемесячного значения в R - PullRequest
Новая
       32

8-дневный MODIS растр до среднемесячного значения в R

0 голосов
/ 24 октября 2018

Я использую набор данных об эвапотранспирации MODIS MOD16A2, который представляет собой сетку 500x500 метров из суммы 8-дневного суммарного испарения.Первый файл заканчивается _001, что означает, что он начинается 1 января и продолжается до 8-го.Следующий файл, _009, затем охватывает период с 8 по 16 января и так далее.

У меня проблемы с переходом от 8-дневных интервалов к среднемесячным значениям.Единственный способ заставить его работать - это найти среднесуточный ET для 8-дневного периода, создать слой на каждый день года, заполнить даты начала среднесуточным ET, перенести последнее наблюдение вперед и разбить на месяцы.оттуда.Проблема в том, что мои растры очень большие, и у меня нет вычислительной мощности или времени, чтобы создать растр глубиной 365 слоев, чтобы объединить его до 12 слоев.

library(raster)

layers <- paste("MODIS", seq(from = 001, to = 365, by = 8), sep = "_")

sampleraster <- brick(nrow = 10, ncol = 10, nl = length(layers))

sampleraster[] <- round(runif(ncell(sampleraster))*50)

names(sampleraster) <- layers

Этот код создает образец растра, похожего на мой, где имена слоев содержат день года, в который начинается период.

Любая помощь в объединении этих 46 слоев с 12 среднемесячными значениями будет принята с благодарностью.

Спасибо

1 Ответ

0 голосов
/ 25 октября 2018

Эта проблема состоит из двух частей.Первый заключается в том, чтобы найти веса, назначаемые каждому слою для вычисления среднемесячных значений.

Часть 1 - Получить весовые коэффициенты

Это не пространственная проблема, и может быть более элегантное решение, но это должно сделать следующее.Обратите внимание, что вам нужно указать год, чтобы иметь возможность справляться с високосными годами.Из вашего описания я понимаю, что последовательности файлов начинаются заново каждый год первого января.

year <- 2017
# number of days in that year (leap year or not?)
ndays <- ifelse(((year %% 100 != 0) & (year %%4 ==0)) | (year %% 400==0), 366 , 365)

# how many layers?
n <- ceiling(ndays/8) 
# day of year for each layer
nn <- rep(1:n, each=8)[1:ndays] 

# day of year for each month
m <- as.integer(format(as.Date(1:ndays, origin=paste0(year-1, "-12-31")), "%m"))

x <- cbind(layer=nn, month=m)

x описывает для каждого дня года, какой слой использовать и какой это месяц.Теперь мы можем для каждого месяца определить, сколько каждого слоя находится в этом месяце (число от 0 до 8 дней). 

weights <- table(x[,1], x[,2])
head(weights)  
 #   1 2 3 4 5 6 7 8 9 10 11 12
 #1 8 0 0 0 0 0 0 0 0  0  0  0
 #2 8 0 0 0 0 0 0 0 0  0  0  0
 #3 8 0 0 0 0 0 0 0 0  0  0  0
 #4 7 1 0 0 0 0 0 0 0  0  0  0
 #5 0 8 0 0 0 0 0 0 0  0  0  0
 #6 0 8 0 0 0 0 0 0 0  0  0  0

Каждый столбец - это месяц, каждая строка - слой.

Часть 2 - Применение весов

Теперь обратимся к примеру растровых данных 

library(raster)
layers <- paste("MODIS", seq(from = 001, to = 365, by = 8), sep = "_")
r <- brick(nrow = 10, ncol = 10, nl = length(layers))
values(r) <- round(runif(ncell(r))*50)
names(r) <- layers

и применим веса для вычисления взвешенных средних для каждого месяца. 

s <- list()
for (i in 1:12) {
    w <- weights[,i]
    x <- r[[which(w > 0)]]
    ww <- w[w > 0] / 8
    s[[i]] <- weighted.mean(x, ww)
}

s <- stack(s)
names(s) <- month.abb
s

#class       : RasterStack 
#dimensions  : 10, 10, 100, 12  (nrow, ncol, ncell, nlayers)
#resolution  : 36, 18  (x, y)
#extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
#names       : Jan, Feb, Mar, Apr, May, Jun, Jul, Aug, Sep, Oct, Nov, Dec 
#min values  :   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0 
#max values  :  50,  50,  50,  50,  50,  50,  50,  50,  50,  50,  50,  50
...