Any workaround to construct temperature distribution over multi-layers raster in R

Question

Here I found a very interesting blog:critical threshold in temperature effects and empirical approach is very interesting, so I want to implement its idea in R. However, I have multi-layer raster data of German' historical daily temperatures (15 years' historical daily mean temperature) in large RasterBrick object. According to the empirical approach that discussed in inspired post, I need to construct the distribution of temperature from my multi-layer raster data.

Update 2: reproducible shapefile:

I aware that downloading shapefile from 3rd party website is not practical, so here I come up reproducible shapefile to give it try:

library(sf)
library(maps)
library(rgeos)
library(mapdata)

germany <- st_as_sf(map("Germany", plot = FALSE, fill = TRUE))
write_sf(germany, "germany.shp")

To easily follow up my post, I created reproducible raster data to work with in R. I also provide Germany' shapefile that taken from eurostat website; here is the shapefile on the fly (I can guarantee the link is quite safe and file is very small to use): eurostat' shapefile and here is handy reproducible raster data:

reproducible data

library(raster)
library(lubridate)
library(tidyverse)

r <- raster(xmn=5.75, xmx= 15, ymn = 47.25, ymx =55,res=c(0.25,0.25))
Deu_crop <- do.call(stack,lapply(1:5479,function(i) setValues(r,round(runif(n = ncell(r),min = -10,max = 25)))))
names(Deu_crop) <- paste0('X',gsub('-','.',ymd('1980.01.01') + days(1:5479)))
shp <- shapefile('eurostat_NUTS3/deu_adm_2006.shp')
e <- raster::extract(Deu_crop,shp)
names(e) <- shp$NUTS_ID

so to test the workflow that presented in inspired post, I need to design several global variables which will be used helper functions that presented in here. But I don't understand how to design some critical global variables that used to accomplish its workflow; It is recommended to define global variable like: w - weather data; tempDat: particular aggregated weather data; Trows: span aggregated grid data; and T: vector of integer temperature (details can be found here: details).

I want to estimate the distribution of temperatures over time from the gridded daily weather data. But I have a difficulty to test the empirical steps that presented in this inspired post because it didn't mention the solution for the case of handling multi-layer raster data, so I don't know how to adopt its fantastic idea on my own in R.

Here is my approach to aggregate multi-layer raster data for each polygons from shapefile (eurostat' shapefile) before using helper funcion in inspired post:

initial attempt to manipulate multi-layers raster:

rasterHelper <- function(ix,e){
  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)) %>% 
    group_by(year) %>% summarise_all(funs(sum)) %>% mutate(NUTS_ID = names(e)[ix])
}
do.call(rbind,lapply(1:length(e),function(ix) rasterHelper(ix)))

but my above attempt is not working; In my attempt, I intend to aggregate temperature raster data for each polygon. The implementation of inspired post is very useful but it is still hard to follow up for handling multi-layer raster data. I assume I should work on each raster layer and construct temperature distribution over time, but really don't have a solid idea how to do in R. Any idea?

Update:

Here is the paper that I got inspiration from it: nonlinear temperature effect ..., but implementing the proposed method is still challenging for me even I followed the workflow that presented in the respective blog: searching critical threshold in temperature effect

Is there anyone point me out how to adopt its empirical approach on multi-layer raster data in R? How can I estimate the distribution of temperature over time? How can I make this happen in R? Any more thoughts? Thanks

Answer 1

I am not quite sure what you want to do.

To set up a smaller example:

library(raster)
lux <- shapefile(system.file("external/lux.shp", package="raster"))
r <- raster(lux)
s <- stack(lapply(1:12, function(i) setValues(r, 1:ncell(r))))

e <- extract(s, lux)

Now you say that you want to aggregate --- that is a bit ambiguous, but perhaps what you want is

x <- lapply(e, function(i) apply(i,2,mean))

equivalent to

y <- extract(s, lux, fun='mean')