Simple linear regression between timeseries raster datasets-CodePudding

I want to run simple linear regression between two (NPP and Tmax) different grid raster datasets. The datasets are annual dataset with similar resolution and time. The first raster dataset is

class      : RasterBrick 
dimensions : 321, 401, 128721, 37  (nrow, ncol, ncell, nlayers)
resolution : 0.0375, 0.0375  (x, y)
extent     : 32.98125, 48.01875, 2.98125, 15.01875  (xmin, xmax, ymin, ymax)
crs        :  proj=longlat  datum=WGS84  no_defs 
source     : Tmax_ANNUAL_RES01.tif 
names      : Tmax_ANNUAL_RES01.1, Tmax_ANNUAL_RES01.2, Tmax_ANNUAL_RES01.3, 
Tmax_ANNUAL_RES01.4, Tmax_ANNUAL_RES01.5, Tmax_ANNUAL_RES01.6, Tmax_ANNUAL_RES01.7, 
Tmax_ANNUAL_RES01.8, Tmax_ANNUAL_RES01.9, Tmax_ANNUAL_RES01.10, Tmax_ANNUAL_RES01.11, 
Tmax_ANNUAL_RES01.12, Tmax_ANNUAL_RES01.13, Tmax_ANNUAL_RES01.14, Tmax_ANNUAL_RES01.15, 
... 
min values :            12.63056,            13.49100,            12.90514,            
11.87769,            12.35822,            12.49712,            12.41137,            
11.83172,            12.46157,             13.08855,             12.18255,             
12.03794,             12.18468,             12.58404,             12.60854, ... 
max values :            40.78825,            39.08575,            41.42342,            
38.99245,            38.76424,            39.27279,            38.86973,            
38.87597,            41.17113,             42.31324,             40.92870,             
42.30246,             42.61915,             43.22092,             41.89821, ...

and The second raster dataset is

class      : RasterBrick 
dimensions : 321, 401, 128721, 37  (nrow, ncol, ncell, nlayers)
resolution : 0.0375, 0.0375  (x, y)
extent     : 32.98125, 48.01875, 2.98125, 15.01875  (xmin, xmax, ymin, ymax)
crs        :  proj=longlat  datum=WGS84  no_defs 
source     : NPP_ANNUAL_RES01.tif 
names      : NPP_ANNUAL_RES01.1, NPP_ANNUAL_RES01.2, NPP_ANNUAL_RES01.3, 
NPP_ANNUAL_RES01.4, NPP_ANNUAL_RES01.5, NPP_ANNUAL_RES01.6, NPP_ANNUAL_RES01.7, 
NPP_ANNUAL_RES01.8, NPP_ANNUAL_RES01.9, NPP_ANNUAL_RES01.10, NPP_ANNUAL_RES01.11, 
NPP_ANNUAL_RES01.12, NPP_ANNUAL_RES01.13, NPP_ANNUAL_RES01.14, NPP_ANNUAL_RES01.15, ... 
min values :                  0,                  0,                  0,                  
0,                  0,                  0,                  0,                  0,                  
0,                   0,                   0,                   0,                   0,                   
0,                   0, ... 
max values :             298.26,             276.47,             281.04,             
274.46,             293.88,             287.74,             289.74,             266.93,             
299.57,              310.09,              304.33,              304.29,              
292.32,              295.15,              283.13, ...

After stacking them the raster datasets is presented as

 stack_ANNUAL <- stack(list.files(pattern="*.tif", full.names=TRUE))

class      : RasterStack 
dimensions : 321, 401, 128721, 74  (nrow, ncol, ncell, nlayers)
resolution : 0.0375, 0.0375  (x, y)
extent     : 32.98125, 48.01875, 2.98125, 15.01875  (xmin, xmax, ymin, ymax)
crs        :  proj=longlat  datum=WGS84  no_defs 
names      : NPP_ANNUAL_RES01.1, NPP_ANNUAL_RES01.2, NPP_ANNUAL_RES01.3, 
NPP_ANNUAL_RES01.4, NPP_ANNUAL_RES01.5, NPP_ANNUAL_RES01.6, NPP_ANNUAL_RES01.7, 
NPP_ANNUAL_RES01.8, NPP_ANNUAL_RES01.9, NPP_ANNUAL_RES01.10, NPP_ANNUAL_RES01.11, 
NPP_ANNUAL_RES01.12, NPP_ANNUAL_RES01.13, NPP_ANNUAL_RES01.14, NPP_ANNUAL_RES01.15, ... 
min values :            0.00000,            0.00000,            0.00000,            
0.00000,            0.00000,            0.00000,            0.00000,            0.00000,            
0.00000,             0.00000,             0.00000,             0.00000,             
0.00000,             0.00000,             0.00000, ... 
max values :          298.26001,          276.47000,          281.04001,          
274.45999,          293.88000,          287.73999,          289.73999,          
266.92999,          299.57001,           310.09000,           304.32999,           
304.29001,           292.32001,           295.14999,           283.13000, ...

The i run the following script to predict NPP from Tmax and want to know the slope

fun=function(x) { if (is.na(x[1])){ NA } else { lm(x[1:37] ~ x[38:74])$coefficients[2] }}
slope <- calc(stack_ANNUAL, fun)

But, the script come up with

Error in (function (classes, fdef, table)  : 
unable to find an inherited method for function ‘writeValues’ for signature 
‘"RasterBrick", "numeric"’

What is the problem?

Thanks!

CodePudding user response：

There is an issue with how you have set up your function for the linear model.

An alternative may be the following:

#since your rasters have the same number of layers we pick any for 'n'
n <- names(NPP) %>% seq_along() %>% max()
#store your independent rasters into separate lists
NPP_list <- c()
TMAX_list <- c()
for( i in 1:n){
  NPP_list[i]<- list(NPP[[i]])
  TMAX_list[i]<- list(TMAX[[i]])
}

#store the rasters into dataframes and remove any na values
NPP_dataframe <- lapply(NPP_list, function(x)data.frame(na.omit(values(x))))
TMAX_dataframe <- lapply(TMAX_list, function(x)data.frame(na.omit(values(x))))

#Map a linear model between the the datasets
mapply(function(x,y)lm(x[[1]] ~ y[[1]], NPP_dataframe, TMAX_dataframe, SIMPLIFY=FALSE)