Let us start by downloading the netCDF file for global daily precipitation for the year 2023 from CHIRPS, accessible through the link: https://data.chc.ucsb.edu/products/CHIRPS-2.0/global_daily/netcdf/p05/chirps-v2.0.2023.days_p05.nc. We can use terra::rast function to import the netCDF to the workspace as a multiband raster. Let us first familiarize ourselves with the dataset.
# Copied path of the rasterdata_path <-"https://data.chc.ucsb.edu/products/CHIRPS-2.0/global_daily/netcdf/p05/chirps-v2.0.2023.days_p05.nc"# Download the raster using download.file. Assign "daily_pcp_2023.nc" name to fileif (file.exists("daily_pcp_2023.nc")==FALSE){download.file(url = data_path, method="curl",destfile ="daily_pcp_2023.nc") }# Plot downloaded filelibrary(terra)pcp=rast("daily_pcp_2023.nc") # Import raster to the environment #~~~ Notice the attributes (esp. nlyr, i.e. number of layers, unit and time)print(pcp)
# Import sample locations from contrasting hydroclimatelibrary(readxl)loc=read_excel("./SampleData-master/location_points.xlsx")print(loc)
# A tibble: 3 × 4
Aridity State Longitude Latitude
<chr> <chr> <dbl> <dbl>
1 Humid Louisiana -92.7 34.3
2 Arid Nevada -116. 38.7
3 Semi-arid Kansas -99.8 38.8
# Value of the lat & lon of the locationslatlon=loc[,3:4] print(latlon)
# A tibble: 3 × 2
Longitude Latitude
<dbl> <dbl>
1 -92.7 34.3
2 -116. 38.7
3 -99.8 38.8
# Plot data for a specific layerworldSHP=terra::vect(spData::world) # Shapefile for CONUSplot(pcp[[100]]) # Same as pcp[[which(time(pcp)=="2023-04-10")]]plot(worldSHP, add=TRUE)points(latlon, pch=19, col="red")
We have previously used terra::extract function to extract cell values from a raster for user-defined locations. However, unlike previous examples, precipitation data is imported from a netCDF, which has 365 layers, one for each day in the year 2023. So, when we use the extract function using point coordinates, 365 values for each location are extracted.
# Import sample locations from contrasting hydroclimatelibrary(readxl)loc=read_excel("./SampleData-master/location_points.xlsx")print(loc)
# A tibble: 3 × 4
Aridity State Longitude Latitude
<chr> <chr> <dbl> <dbl>
1 Humid Louisiana -92.7 34.3
2 Arid Nevada -116. 38.7
3 Semi-arid Kansas -99.8 38.8
# Value of the lat & lon of the locationslatlon=loc[,3:4] print(latlon)
# A tibble: 3 × 2
Longitude Latitude
<dbl> <dbl>
1 -92.7 34.3
2 -116. 38.7
3 -99.8 38.8
# Extract time series using "terra::extract"loc_pcp=terra::extract(pcp, latlon, #2-column matrix or data.frame with lat-longmethod='bilinear') # Use bilinear interpolation (or ngb) option# View data sampleloc_pcp[,1:8] # View(loc_pcp)
# Plot hyetograph for the location in Louisianalibrary(ggplot2)pcp_df1=data.frame(time=time(pcp),pcp=as.numeric(loc_pcp[1,-c(1)])) # Select first row, exclude the first column# ggplotggplot(pcp_df1,aes(x=time,y=pcp)) +geom_bar(stat ='identity')+theme_bw()+ylab("Precipitation [mm/day]")+xlab("Time [Days]")
# Export the extracted data as CSVwrite.csv(loc, "extracted_pcp2023.csv")
8.3 Cell–wise Operation on All Layers
The terra::app() function applies a function to each cell of a raster and is used to summarize (e.g., calculating the sum) the values of multiple layers into one layer. We will use the NDVI netCDF generated in Sec 7.1.2 for the year 2016 to calculate global cellwise statistics.
#Let's look at the help section for app()# ?terra::app# Import NDVI layers from the NDVI generated in 7.1.2ndviStack =rast("NDVI.nc")ndviStack
# Subset raster stack/brick (notice the double [[]] bracket and similarity to lists)ndviStack_sub=ndviStack[[c(1,3,5,10,12)]] #Select 1st, 3rd, 5th, 10th and 12th layers# Calculate mean of each grid cell across all layersmean_ras =app(ndviStack_sub, fun=mean, na.rm = T)# Calculate sum of each grid cell across all layerssum_ras =app(ndviStack_sub, fun=sum, na.rm = T)#~~ A user-defined function for mean, variance and skewnessmy_fun =function(x){ meanVal=mean(x, na.rm=TRUE) # Mean varVal=var(x, na.rm=TRUE) # Variance skewVal=moments::skewness(x, na.rm=TRUE) # Skewness output=c(meanVal,varVal,skewVal) # Combine all statisticsnames(output)=c("Mean", "Var","Skew") # Rename output variablesreturn(output) # Return output} # Apply user function to each cell across all layersstat_ras =app(ndviStack_sub, fun=my_fun)# Plot statisticscoastlines =vect("./SampleData-master/ne_10m_coastline/ne_10m_coastline.shp")library(cetcolor)colpal = cetcolor::cet_pal(20, name ="r2") plot(stat_ras, col =rev(colpal), # Rev argument reverses the Palette (from Bu-> Rd to Rd -> Bu)asp =NA, # Aspect ratio: NA, fill to plot spacenc =2, # Number of columns to arrange plotsfun =function(){plot(coastlines, add=TRUE)} # Add coastline boundary)
8.4 Cell–wise Operation on Layer Groups
terra::tapp() is an extension of app(), allowing us to select a subset of layers for which we want to perform a certain operation. Let’s have the first two layers as group 1 and the next three as group 2. Function will be applied to each group separately and 2 layers of output will be generated.
#The layers are combined based on indexing.stat_ras = terra::tapp(ndviStack_sub,index=c(1,1,2,2,2),fun= mean)# Try other functions: "sum", "mean", "median", "modal", "which", "which.min", "which.max", "min", "max", "prod", "any", "all", "sd", "first".names(stat_ras) =c("Mean_of_rasters_1_to_2", "Mean_of_rasters_3_to_5")# Two layers are formed, one for each group of indices# Lets plot the two output rastersplot(stat_ras, col =rev(colpal), # Rev argument reverses the Palette (from Bu-> Rd to Rd -> Bu)asp =1, # Aspect ratio: 1, i.e. X==Ync =1, # Number of columns to arrange plotsfun =function(){plot(coastlines, add=TRUE)} # Add coastline boundary)
8.5 Layers as Function Arguments
The terra::lapp() function allows to apply a function to each cell using layers as arguments.
#Let's look at the help section for app()# ?terra::lapp#User defined function for finding differencediff_fun =function(a, b){ return(a-b) }diff_rast =lapp(ndviStack_sub[[c(4, 2)]], fun = diff_fun)#Plot NDVI differenceplot(diff_rast, col =rev(colpal), # Rev argument reverses the Palette (from Bu-> Rd to Rd -> Bu)asp =1, # Aspect ratio: 1, i.e. X==Yfun =function(){plot(coastlines, add=TRUE)} # Add coastline boundary)
We will use a time series of global gridded NDVI rasters to fit a separate linear regression model at every land pixel. By the end of this lesson, you should be able to:
Read multiple raster files as a time-ordered SpatRaster using the terra package.
Understand how a raster stack represents a time series at each grid cell.
Create lagged NDVI predictors.
Fit a multiple linear regression model for one pixel.
Apply the model independently to every raster cell with terra::app().
Handle non-land pixels and model failures using NA values and tryCatch().
Interpret and map the resulting regression coefficients.
8.6 Cell-wise regression and prediction
Let us consider a grid cell with NDVI observations through time:
\[
NDVI_1, NDVI_2, NDVI_3, \ldots, NDVI_T.
\]
We model current NDVI as a function of its previous three layers:
Here, \(\beta_0\) is the intercept; \(\beta_1\), \(\beta_2\), and \(\beta_3\) describe the relationships with NDVI one, two, and three layers earlier; and \(\epsilon_t\) represents unexplained variation. This is an autoregressive-style regression because NDVI values from previous layers are used to explain current NDVI. The model is implemented separately for every grid cell.
Important
A lag of one means one raster layer, not necessarily one day or one month. Its physical meaning depends on the temporal spacing of the input files.
8.6.1 Prepare the multilayer NDVI raster
We will use the terra package to import the raster layers, inspect the resulting SpatRaster, and apply a function across the time series of each grid cell.
library(terra)
First, locate the NDVI files and combine them into a multilayer raster.
ndvi_directory <-paste0("./SampleData-master/raster_files/NDVI")rasNames <-list.files(path = ndvi_directory,pattern ="\\.tif$",full.names = T)# The regression requires the rasters to be in chronological order.rasNames <-sort(rasNames)ndviStack <-rast(rasNames)ndviStack
list.files() searches the directory for files ending in .tif. Setting full.names = T returns the complete file paths required by rast(). The sort() function places filenames in alphanumeric order. This produces chronological order only when dates in the filenames use a sortable format, such as NDVI_2020_01.tif and NDVI_2020_02.tif. Inspect the filenames before fitting the model.
plot(ndviStack[[1]], main ="NDVI at the first time step")
For any grid cell, its values across all layers form an NDVI time series. Let us extract one example cell to see the type of data passed to the cellwise function.
# Choose a cell number known to be on land for the sample data.example_cell <-1000x_example <-as.numeric(ndviStack[example_cell])x_example
When app() evaluates a cell, it passes a vector like x_example to the function. Non-land cells contain only NA values. Land cells may also contain occasional NA values because of clouds, missing observations, or quality screening.
8.6.2 Construct the lagged regression model
Let us first construct lagged predictors for a short example time series.
current_NDVI lag_1 lag_2 lag_3
[1,] 0.20 NA NA NA
[2,] 0.25 0.20 NA NA
[3,] 0.31 0.25 0.20 NA
[4,] 0.29 0.31 0.25 0.20
[5,] 0.35 0.29 0.31 0.25
[6,] 0.39 0.35 0.29 0.31
The first row cannot have a one-layer lag, so lag_1 begins with NA. Similarly, lag_2 begins with two missing values and lag_3 begins with three. Therefore, the first three rows cannot be used to fit the regression.
We can now place these steps inside a function that receives the NDVI time series for one cell and returns the four model coefficients.
NDVI_pred_mod <-function(x) {# Convert the cell values to a numeric vector. x <-as.numeric(x)tryCatch({# Skip cells with no NDVI observations, such as non-land pixels.if (!all(is.na(x))) {# Construct one-, two-, and three-layer lagged predictors. x1 <-c(NA, x[-length(x)]) x2 <-c(NA, x1[-length(x1)]) x3 <-c(NA, x2[-length(x2)])# Column 1 is current NDVI; columns 2-4 are lagged NDVI. lagInput <-cbind(x, x1, x2, x3)# Fit current NDVI using the previous three layers. regCoeff <-lm(lagInput[, 1] ~ lagInput[, 2:4])# Extract the intercept and slope coefficients. interc <- regCoeff$coefficients[1] a1 <- regCoeff$coefficients[2] a2 <- regCoeff$coefficients[3] a3 <- regCoeff$coefficients[4]return(c(interc, a1, a2, a3)) } else {return(rep(NA, 4)) } }, error =function(e) {return(rep(NA, 4)) })}
The function must always return four values: the intercept and three lag coefficients. Non-land cells and failed models also return four values, but each value is NA. This keeps the output structure consistent across the raster.
tryCatch() prevents one problematic cell from stopping the entire operation. A model may fail because it has too few complete observations, no variation in valid NDVI values, or strong dependence among the lagged predictors.
Important
tryCatch() handles errors, but it can also hide their causes. While developing the function, test it on selected cells without tryCatch() so that you can examine any model errors.
Before applying the function globally, test it on a valid example and an all-NA cell.
set.seed(42)test_series <-0.5+0.2*sin(seq(0, 4* pi, length.out =30)) +rnorm(30, mean =0, sd =0.02)plot(test_series, type="l", main="sample auto-correlated time series")
The first call should return four coefficients. The second should return four NA values.
8.6.3 Apply and interpret the cellwise model
We will now run the function over the complete global raster. During rendering, app() passes the NDVI time series of every grid cell to NDVI_pred_mod() and creates four global coefficient layers.
terra::app() gathers the values across all layers at one grid cell, supplies that vector to the function, and writes the four returned coefficients to four output raster layers. Because this code chunk is evaluated, the global coefficient rasters are created before the plotting chunk runs. The output retains the geometry of the input raster.
# Plot all four coefficient layers in the rendered output.plot_colors =colorRampPalette(c("saddlebrown", "gold", "lightyellow", "lightblue", "darkblue"))(100)# Arrange the four coefficient maps in a 2 x 2 layoutpar(mfrow =c(2, 2))plot(coeffVal[["intercept"]],main ="Intercept",range =c(0, 1),col = plot_colors)plot(coeffVal[["lag_1_coefficient"]],main ="Lag-1 coefficient",range =c(0, 3),col = plot_colors)plot(coeffVal[["lag_2_coefficient"]],main ="Lag-2 coefficient",range =c(-2, 2),col = plot_colors)plot(coeffVal[["lag_3_coefficient"]],main ="Lag-3 coefficient",range =c(-1, 1),col = plot_colors)
# Restore the default graphics layoutpar(mfrow =c(1, 1))
A positive lag coefficient indicates that, after accounting for the other lags, higher NDVI in an earlier layer is associated with higher current NDVI. A negative coefficient indicates an inverse partial relationship. Larger absolute coefficients represent stronger fitted relationships on the scale of the model. The intercept is the predicted NDVI when all three lagged NDVI values equal zero, although that condition may not be ecologically meaningful everywhere.
A coefficient map is not a prediction map. It shows how fitted model parameters vary spatially. For one time series, the next NDVI value can be estimated from the three most recent observations:
Why do the first three rows of lagInput contain missing values?
What does a lag of one represent for this NDVI dataset?
Why must the raster files be arranged chronologically?
Why does the function return four values for a non-land pixel?
What happens if a land pixel has too few complete lagged observations?
Does a positive lag-1 coefficient prove that earlier NDVI causes current NDVI? Why or why not?
How would you modify the model to use only two lags?
How would you produce a raster of predicted NDVI for the next time step?
Why might neighboring pixels have different coefficients?
What diagnostics would you inspect before using the model for forecasting?
The key elements of this workflow are the temporal ordering of the raster layers, the construction of lagged predictors, consistent handling of missing values, and the use of terra::app() to repeat the same analysis independently at every grid cell.
