Concatenation (placing the objects one after the other) of rasters follow the same syntax as arrays i.e. by using c(obj1, obj2). Concatenation of rasters can be used for, but not limited to, creating a raster time series to evaluate time-dependent variability in the spatial data. Here we will use the global NDVI rasters from MODIS (resampled to a coarse resolution for computational ease) to create a raster time series in the form of a multilayer rast object, and will demonstrate several spatial and arithmetic operations on the multilayer rasters.
7.2.1 Create Raster Time Series
The folder ./SampleData-master/raster_files/NDVI/ contains 23 rasters of global NDVI for the year 2016, with a retrieval frequency of 16 days. Let us import four sample rasters and concatenate them to generate a multilayer raster.
#~~~ Create and plot NDVI SpatRasterlibrary(terra)library(cetcolor)# Location of the NDVI raster filesndvi_path="./SampleData-master/raster_files/NDVI/"# List of all NDVI rastersras_path=list.files(ndvi_path, # The folder which contains the files of interestpattern='*.tif', # Select files with this extentionfull.names=TRUE) # Print full file path? Yes (TRUE) or no (FALSE)?print(ras_path) # The folder contains 23 rasters, in increasing order of time
# User the file path to import rasterr1=rast(ras_path[1])r2=rast(ras_path[2])r3=rast(ras_path[3])r4=rast(ras_path[4])# Concatenate the rasters i.e. arrange the rasters one after the otherrast_ts=c(r1,r2, r3, r4) nlyr(rast_ts) # Number of layers in the multilayer raster
# Plot concatenated rasterscolpal = cetcolor::cet_pal(20, name ="r2") plot(rast_ts, col=rev(colpal))
7.2.2 Programmatic Concatenation of rast Objects
Concatenation method shown above can be good for a few rasters. However, in the case of a long time series, this method could be cumbersome. In such cases, we can write a small script to concatenate the rasters in to a multilayer rast object. We will first explore the use of loops- which will concatenate each rast layer sequentially. This method may take longer, but helps to avoid memory bottleneck when a very large number of rasters are simultaneously processed (the threshold may depend on the system memory and processor). In contrast, lapply and map vectorize the operation of creating the rast object for concatenation, and maybe faster for a moderate number of files (again, depends on the system configuration).
#~~~~~~~~~~~~~~~~~~# Method 1: Using for loop to create raster layers from the raster location#~~ Suitable for large number of rasters and avoiding memory bottleneckras_stack=rast()for (nRun in1:length(ras_path)){ ras_stack=c(ras_stack,rast(ras_path[[nRun]]))}# Check dimension of data cubedim(ras_stack) #[x: y: z]- 23 raster layers with 456 x 964 cells
[1] 456 964 23
#~~~~~~~~~~~~~~~~~~# Method 2: Use lapply to create raster layer list from the raster file pathsras_list =lapply(paste(ras_path, sep =""), rast)# This a list of 23 raster objects stored as individual elements.# Convert raster layer lists to Spatraster ras_stack =rast(ras_list) # Stacking all rasters as a SpatRaster # This a multi-layer (23 layers in this case) SpatRaster Object.#~~~~~~~~~~~~~~~~~~# Method 3: maping the rast function on the raster file pathslibrary(purrr) ras_list = purrr::map(ras_path, ~rast(.x)) # Import the rasters into a listras_stack =rast(ras_list) # Convert a list of rasters to rast objectdim(ras_stack) #[x: y: z]- 23 raster layers with 456 x 964 cells
[1] 456 964 23
Computing Counsel: Exporting concatenated rasters as a netCDF is a convenient time-saving measure. We can use terra::writeCDF to export SpatRaster object as netCDF. This allows the multilayer raster to be imported directly to the workspace and avoid repeating the steps in the preceeding section.
# Export concatenated rasters as a netCDFterra::writeCDF(ras_stack, # SpatRasterfile.path("NDVI.nc"), # Output filenameoverwrite=TRUE,varname="NDVI",unit="[-]",longname="Global MODIS NDVI, 2016, 16_day, 36KM",zname='time')
7.3 Spatial Operations on Multilayer Raster
7.3.1 Subset Multilayer Raster
Subsetting multilayer rasters (i.e. creating of smaller subset of layers) follows the same syntax we would use to subset a list. Notice the double square bracket i.e. [[ ]] we use to subset raster layers. When generating a spatial plot of multilayer rasters using terra::plot, to add vector boundaries to all plot layers, we need to define a custom function to add the object to each layer.
# Subset raster stack/brick (notice the double [[]] bracket and similarity to lists)sub_ras_stack=ras_stack[[c(1,3,5,10,12)]] #Select 1st, 3rd, 5th, 10th and 12th layers#~~ Plot first 4 elements of NDVI SpatRastercoastlines =vect("./SampleData-master/ne_10m_coastline/ne_10m_coastline.shp")# Function to add shapefile to the plotaddCoastlines=function(){plot(coastlines, add=TRUE) # Add Coastline vector to existing plot } # Notice how a subset of first 4 rasters is created terra::plot(sub_ras_stack[[1:4]], # Select raster layers to plotcol=colpal, # Specify colormapasp=NA, # Aspect ratio= NAfun=addCoastlines) # Add coastline to each layer
7.3.2 Crop and Mask
Cropping and masking a multilater raster follows the same format as that of a single raster we learned in Sec 5.3.
# Names of the states to mask NDVIstates =c('Oklahoma','Texas','New Mexico', 'Louisiana', 'Arkansas')# Subset the selected states from CONUS simple feature (sf) objectlibrary(sf)library(spData)us_shp=spData::us_states# Subset sf for selected states using the NAME featureselectState = us_shp[(us_shp$NAME %in% states),]# Convert selectState sf to vect objectstateVect=vect(selectState)# Reproject the state vector to match the CRS of the NDVI rasterstateVect_proj=terra::project(stateVect,crs(ras_stack))# Raster operationndvi_crp =crop(ras_stack, ext(stateVect_proj)) # Crop rasterndvi_msk = terra::mask(ndvi_crp,stateVect_proj) # Maskdim(ndvi_msk) #Notice the number of layers in the masked rast
[1] 30 53 23
# Plot raster and shapefileplot(ndvi_msk[[1:4]], col=colpal, fun=function(){plot(stateVect_proj, add=TRUE)} # Add state borders )
7.3.3 Spatial Extraction
Similar to single-band raster operation shown in Sec. 5.5.1.2, this operation can be executed using terra::extract. However, the output will contain values from all layers in the multiband raster.
# Extract cell values within each feature for all layersndvi_state_mean = terra::extract(ndvi_msk, # Data cube stateVect_proj, # Shapefile for feature reference na.rm=T)# Extract 'mean' statistic of cell values within each feature for all layersndvi_state_mean = terra::extract(ndvi_msk, # Multilayer rast object stateVect_proj, # Vect object of statesfun=mean, # Summary statistics (mean/sum/min/max)na.rm=T ) # Remove NA values? # Try: View(ndvi_state_mean)#~~~ Each row corresponds to one state. Columns store layer-wise mean NDVI per statehead(ndvi_state_mean)
7.4 Layer–wise Statistical/Arithmetic Operations on Multilayer Raster
7.4.1 Layer–wise Statistics
Similar to what we learned in Sec. 4.2 Raster Statistics, we can use the global function to obtain summary statistics of the multilayer raster. Here the operation will yield the desired statistic for each constituent layer.
# Layer-wise cell-statistics # Layer-wise Meanglobal(sub_ras_stack, mean, na.rm= T) # Try modal, median, min etc.
# User-defined statistics by defining own functionquant_fun =function(x) {quantile(x, probs =c(0.25, 0.75), na.rm=TRUE)} global(sub_ras_stack, quant_fun) # 25th, and 75th percentile of each layer
Again, this operation is similar to the operations on a single-band rasters, except the operation is carried out for each constituent layer simultaneously. Let’s explore some examples.
# Arithmetic operations on SpatRaster are same as listsadd = ndvi_msk+10# Add a number to raster layersmult = ndvi_msk*5# Multiply a number to raster layerssubset_mult = ndvi_msk[[1:3]]*10# Multiply a number to a subset of raster layerslog_ra =log(ndvi_msk) # Layer-wise Log-transformation# Data filtering based on cell-value#~~~ Assign NA to any value less than 0.5#~~~ Try ?terra::clamp for help and other options/usesfilter_rast = terra::clamp(ndvi_msk, lower=0.5, values=FALSE)# Let's plot the filtered rastersplot(filter_rast[[1:4]], asp=NA, # Aspect ratio: NA, fill to plot spacecol=colpal, # Assign colormap to the plotnc=2, # Number of columns to arrange plots, range=c(0,0.9), # Set custom z range for the plotsfun=function(){plot(stateVect_proj, add=TRUE)} # Add state borders)
# Summary statistics of layers using boxplots#~~~ Type ?boxplot in console for details on the selected options below#~~~ Based on the boxplot, can we be sure that clamping worked? dateStr =substr(names(ndvi_msk),13,23) # Substring of dates for X-axis labelsboxplot(filter_rast, horizontal=F, col ="royalblue",frame=F, cex.axis=1, notch=T, col.axis ="black",whisklty=1,whiskcol="black", whisklwd=1,staplelwd =1, staplecol ="royalblue",medcol ="red",medlwd=1.2, col.ticks ="black", outlier=F, names=dateStr) axis(1, tck=0.05, col.axis ="transparent", lwd=1)
# Layerwise normalization to [0,1]#~~~ Calculate minumum and maximum cell values for each layermin_val =global(ndvi_msk, min, na.rm= T) # Layer-wise minimamax_val =global(ndvi_msk, max, na.rm= T) # Layer-wise maxima#~~~ Normalize layers using respective minima and maximanorm_stk=(ndvi_msk-min_val$min)/(max_val$max-min_val$min) # Plot the normalized rasters. Notice the raster values range between 0 and 1plot(norm_stk[[1:4]],col=colpal)
