Вот решение (включая воспроизводимый пример):
library(raster)
library(lubridate)
library(tidyverse)
# creating some fake temperature data which matches your rasterstack
# create template raster
r <- raster(xmn=5.75, xmx= 15, ymn = 47.25, ymx =55,res=c(0.25,0.25))
# add fake temperature values
Deu_crop <- do.call(stack,lapply(1:5479,function(i) setValues(r,round(runif(n = ncell(r),min = -10,max = 25)))))
# add layer names
names(Deu_crop) <- paste0('X',gsub('-','.',ymd('1980.01.01') + days(1:5479)))
# check rasterstack
Deu_crop
# output
#
# class : RasterStack
# dimensions : 31, 37, 1147, 5479 (nrow, ncol, ncell, nlayers)
# resolution : 0.25, 0.25 (x, y)
# extent : 5.75, 15, 47.25, 55 (xmin, xmax, ymin, ymax)
# coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
# names : X1980.01.02, X1980.01.03, X1980.01.04, X1980.01.05, X1980.01.06, X1980.01.07, ...
# min values : -10, -10, -10, -10, -10, -10, ...
# max values : 25, 25, 25, 25, 25, 25, ...
Так что Deu_crop должно быть сопоставимо с вашими данными с точки зрения структуры, конечно со случайными значениями температуры.
Шейп-файл нелегко воспроизвести, поэтому я скачал ваш и использовал его.Как я уже упоминал, некоторые полигоны немного малы для извлечения.
Самый быстрый способ сделать это - растеризовать шейп-файл в соответствии с растром данных, но некоторые полигоны не будут преобразованы, а другие, вероятно, в неправильную ячейку ... так что в этом случае это может бытьЛучше использовать raster::extract напрямую с шейп-файлом, даже если он немного медленный.Но если вам нужно сделать это всего пару раз, это сносно - выпейте кофе тем временем.
shp <- shapefile('eurostat_NUTS3_29-May-18/deu_adm_2006.shp')
# coffee time
e <- extract(Deu_crop,shp)
# add NUTS_ID as names to list
names(e) <- shp$NUTS_ID
Чтобы рассчитать количество дней в году на одну корзину, я создаю функцию, которая использует tidiverse функциональность и использование lapply для перебора всего списка извлечений (один элемент списка соответствует одному многоугольнику):
# define bins
bins <- seq(-10,25,length.out = 5)
myfun <- function(ix){
gather(data.frame(e[[ix]],stringsAsFactors = F),'colname','temp') %>%
group_by(colname) %>% summarise(temp = mean(temp)) %>% ungroup() %>% # spatial mean
mutate(year = sub('X(\\d{4}).+','\\1',colname)) %>% # get years
select(- colname) %>% # drop colname column
mutate(bin1= (temp <= bins[1]) * 1) %>% # bin1
mutate(bin2= (temp > bins[1] & temp <= bins[2]) * 1) %>% # bin2
mutate(bin3= (temp > bins[2] & temp <= bins[3]) * 1) %>% # bin3
mutate(bin4= (temp > bins[3] & temp <= bins[4]) * 1) %>% # bin4
mutate(bin5= (temp > bins[4] & temp <= bins[5]) * 1) %>% # bin5
mutate(bin6= (temp > bins[5]) * 1) %>% select(-temp) %>% # bin6
group_by(year) %>% summarise_all(funs(sum)) %>% mutate(NUTS_ID = names(e)[ix]) # drop year, calculate occurences and add NUTS_ID
}
# create single dataframe
result <- do.call(rbind,lapply(1:length(e),function(ix) myfun(ix)))
Быстрый просмотр переменной result:
result
# output:
#
# # A tibble: 6,864 x 8
# year bin1 bin2 bin3 bin4 bin5 bin6 NUTS_ID
# <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <chr>
# 1 1980 12 85 91 92 85 0 DEA54
# 2 1981 3 64 99 113 86 0 DEA54
# 3 1982 3 80 113 86 83 0 DEA54
# 4 1983 6 84 90 85 100 0 DEA54
# 5 1984 8 90 92 86 90 0 DEA54
# 6 1985 5 86 85 95 94 0 DEA54
# 7 1986 6 74 97 108 80 0 DEA54
# 8 1987 4 82 99 94 86 0 DEA54
# 9 1988 3 89 87 91 96 0 DEA54
#10 1989 8 103 92 73 89 0 DEA54
# # ... with 6,854 more rows
Обновление:
Для обработки бинов я сначала вычисляю бины по минимуму и максимуму всех данных, а затем использую новую функцию createBins, чтобы добавить их к извлечению каждого многоугольника.Это заменит myfun часть моего оригинального решения.
# new function
createBins <- function(df,bins_mat){
for (i in 1:nrow(bins_mat)){
bin <- sprintf('Bin%s;%s;%s',bins_mat[i,1],bins_mat[i,2],bins_mat[i,3])
if (i ==1) df <- df %>% mutate(!!bin := (temp >= bins_mat[i,2] & temp <= bins_mat[i,3])*1)
else df <- df %>% mutate(!!bin := (temp > bins_mat[i,2] & temp <= bins_mat[i,3])*1)
}
return(df)
}
# new version of myfun
myfun2 <- function(ix,bins_mat){
gather(data.frame(e[[ix]],stringsAsFactors = F),'colname','temp') %>%
group_by(colname) %>% summarise(temp = mean(temp)) %>% ungroup() %>% # spatial mean
mutate(year = sub('X(\\d{4}).+','\\1',colname)) %>% # get years
select(- colname) %>% # drop colname column
createBins(.,bins_mat) %>% select(-temp) %>%
group_by(year) %>% summarise_all(funs(sum)) %>% mutate(NUTS_ID = names(e)[ix])
}
# 11 values to create 10 interval bins
bins <- seq(min(cellStats(Deu_crop,'min')),min(cellStats(Deu_crop,'max')),length.out = 11)
# create a bin matrix (number, bin_minimum, bin_maximum) for later function
bins_mat <- cbind(1:10,bins[1:10],bins[2:11])
# create new result
result <- do.call(rbind,lapply(1:length(e),function(ix) myfun2(ix,binsmat)))