Any efficient workaround to reshape and aggregate multi-layers raster grid in R?

Question

I have historical daily minimum/maximum temperature observations data in two big netCDF file, so I used raster::stack to import those in RasterStack object. However, here I want to reshape and aggregate two different multi-layers raster grid into one plain text ASCII file (want to replicate of this).

reproducible data

Here is the reproducible example for historical daily minimum and daily maximum temperature observation in multi-layers raster grid:

T_max <- raster(xmn=5.75, xmx= 15, ymn = 47.25, ymx =55,res=c(0.25,0.25))
Deu_df_max <- do.call(stack,lapply(1:5479,function(i) setValues(T_max,round(runif(n = ncell(T_max),min = 2,max = 35)))))
names(Deu_df_max ) <- paste0('X',gsub('-','.',ymd('1980.01.01') + days(1:5479)))

T_min <- raster(xmn=5.75, xmx= 15, ymn = 47.25, ymx =55,res=c(0.25,0.25))
Deu_df_min <- do.call(stack,lapply(1:5479,function(i) setValues(T_min ,round(runif(n = ncell(T_min ),min = -5,max = 24)))))
names(Deu_df_min ) <- paste0('X',gsub('-','.',ymd('1980.01.01') + days(1:5479)))

Here is what I tried:

To do so, I extracted German' shapefile from Eurostat website, here you could see minimal shapefile on the fly (it is perfectly safe, and if you don't want to download, I will provide reproducible shapefile to use). Please take a quick look original minimal shapefile. And I did raster extraction as follow:

library(raster)
shp <- shapefile('eurostat_NUTS3_29-May-18/deu_adm_2006.shp')
deu_extr_mx <- raster::extract(Deu_df_max ,shp)
deu_extr_mn <- raster::extract(Deu_df_min ,shp)
names(deu_extr_mx ) <- shp$NUTS_ID
names(deu_extr_mx ) <- shp$NUTS_NAME
names(deu_extr_mn ) <- shp$NUTS_ID
names(deu_extr_mn ) <- shp$NUTS_NAME

based on the output of deu_extr_mx and deu_extr_mn, I intend to aggregate grid temperature data for each polygon in each year respectively.

Basically, I want to reshape and aggregate T_max and T_min raster grid for each polygon and do merge in plain text ASCII format. I don't have any solid idea how to make the aggregation and reshape data for multi-layers raster grid in R.

Can anyone give me any possible idea or programmatic approach to make this happen in R? How can I reshape and aggregate the raster grid data for the respective polygon and put the merging results in a simple plain text? Any idea to get this done in R more efficiently?

desired output:

I want to aggregate grid temperature observations for each polygon (minimal shapefile on the fly) in each year respectively and put the aggregation result in plain text ASCII format. Is that doable?

Here is the skeleton of my desired aggregation output in plain text ASCII format (I fill it with the random number):

 date   year month day NUTS_ID  NUTS_NAME        lat         long    tmax        tmin
1 1980-1-1 1980    1  1  DE111      Hamm            38.5        -122.5  12.5000     4.1100
2 1980-1-2 1980    1  2  DE112      Ostvorpommern   38.5        -122.5  12.5445     4.2894
3 1980-1-3 1980    1  3  DE113      Aurich          38.5        -122.5  12.5878     4.4574
4 1980-1-4 1980    1  4  DE114      Ludwigslust     38.5        -122.5  12.6298     4.6144
5 1980-1-5 1980    1  5  DE115      Wilhelmshaven   38.5        -122.5  12.6706     4.7604
6 1980-1-6 1980    1  6  DE116      Mecklenburg     38.5        -122.5  12.7103     4.8959

Is that doable to reshape, merge and aggregate temperature grid in multi-layers raster grid for each polygons in each year in R? Is there any efficient way to make this happen? Any idea?

Answer 1

Overall the way you are wanting to approach your analysis is not a good idea and is loosing the overlay and memory-safe advantages of the raster class. It is not uncommon, for those not familiar with the raster format, to go down the road of a flat file (table) format. With "smallish" rasters, this is doable but when you get into high-dimension data the dimensionality can explode quite quickly. For instance, you have a climate raster of rows = 1000 & columns = 1000 (n=1000,000), which is seemingly a tractable dataset. However, if you add a t dimension of say, days across 30 years (n=10,950), you suddenly end up with n = 10,950,000,000. Just food for thought.

You can create the format you are after by coercing your raster stack data into an sp SpatialPixelsDataFrame object (or SpatialPointsDataFrame using raster::rasterToPoints). We can work though an example starting with creating some dummy data for tmin and tmax.

library(raster)

f <- raster(nrows=50, ncols=50, xmn=0, xmx=10)
tmax <- stack()
  for(i in 1:12) {
    x <- f
    x[] <- runif(ncell(x),20,90)
    tmax <- addLayer(tmax, x)
}   
tmin <- stack()
  for(i in 1:12) {
    x <- f
    x[] <- runif(ncell(x),20,90)
    tmin <- addLayer(tmin, x)
}       
( d <- seq(as.Date("2000/1/1"), as.Date("2000/12/31"), "months") )

names(tmin) <- paste("tmin",d,sep="_")
names(tmax) <- paste("tmax",d,sep="_")

Now we can coerce the data into a SpatialPixlesDataFrame object where the @data slot will contain a data.frame object with the attributes.

temp <- as(tmin, "SpatialPixelsDataFrame")
tmax.sp <- as(tmax, "SpatialPixelsDataFrame")

Since the data represent the same arrays, you can just directly combine the two data.

temp@data <- data.frame(temp@data, tmax.sp@data) 
   head(temp@data)

You can now use something like sp::over to assign attributes from polygon data. You can pull this data into a data.frame (or other class) object by assigning the @data slot (and even add coordinates).

temp <- data.frame(coordinates(temp), temp@data)