vector-from-raster

library(msdrought)

In general, the msdrought Package works best when the data from a specific point are a raster is pulled, converted into an xts object, and then analyzed. This vignette will show an alternative approach, when an entire raster is analyzed for its duration and intensity values, then a single specific point and its corresponding data is extracted to a vector. The process for both duration and intensity (as well as any other given parameter for the msdStats function) is identical.

# Load the necessary packages and the included data, for the purpose of demonstration
library(terra)
library(tidyr)
library(lubridate)
library(stringr)
library(ggplot2)
library(xts)

# This loads the data included in the package, but you would attach your own
data <- system.file("extdata", "prcp_cropped.tif", package = "msdrought") 
infile <- terra::rast(data)

Now we have our raster data loaded in as “infile”. Using the terra::app command, we can apply a function to each cell of the data. Begin by finding the range of dates for the raster data, as these values are needed for the package to understand the statistics. Data processing begins by filtering the data using a bartlett noise filter, then apply the msdStats function to the data.

# Find the key dates related to the MSD
# msdDates = msdDates(x, firstStartDate, firstEndDate, secondStartDate, secondEndDate)
formattedDates <- as.Date(terra::time(infile))
critDates <- msdrought::msdDates(formattedDates)

# Use the terra::app function to apply the Bartlett noise filter (msdFilter) to the raster
# msdFilter = msdFilter(x, window)
filtered <- suppressWarnings(terra::app(infile, msdFilter, window = 31, quantity = 2))
terra::time(filtered) <- formattedDates

# Use the terra::app function to apply the Bartlett noise filter (msdFilter) to the raster
# msdStats = msdStats(x, dates, fcn)
intensityPlots <- suppressWarnings(terra::app(filtered, msdStats, critDates, fcn = "intensity"))
durationPlots <- suppressWarnings(terra::app(filtered, msdStats, critDates, fcn = "duration"))
# Note: suppressWarnings hides irrelevant messages that result from the msdFilter output

The “intensityPlots” and “durationPlots” rasters contain the intensity and duration values for every geographical location for both of the years of data contained in the original dataset. Now that we have all the intensity and duration values, we will extract a single point’s data from the new rasters.

# Set up the desired location's longitude and latitude values
lon <- -86.2621555581 # Longitude of the spatial point we're interested in analyzing
lat <- 13.3816217871 # Latitude of the spatial point we're interested in analyzing
lonLat <- data.frame(lon = lon, lat = lat)

# Set up precipitation data by extracting the data located at our longitude and lattitude
location <- vect(lonLat, crs = "+proj=longlat +datum=WGS84")
locationIntensity <- terra::extract(intensityPlots, location, method = "bilinear") |>
  subset(select = -ID) |>
  t()

locationDuration <- terra::extract(durationPlots, location, method = "bilinear") |>
  subset(select = -ID) |>
  t()

Now we have converted the “intensityPlots” and “durationPlots” rasters to numeric objects called “locationIntensity” and “locationDuration”, which represent the intensity and duration values for the provided coordinates across both years of data. If a different location’s data were desired, all we would need to do is swap out the lon and lat values we provided and run the previous chunk of code again. For ease of comparison to the “Sample-Walkthrough” vignette, we’ll combine the two objects together. If we had repeated the process for all 5 msdStats parameters (“duration”, “intensity”, “firstMaxValue”, “secondMaxValue”, or “min”), we could achieve the same output as the “Sample-Walkthrough” vignette.

# Pull the years from the data, for the purpose of organizing the data
allYears <- terra::time(infile)
firstYear <- lubridate::year(allYears[1])
lastYear <- lubridate::year(allYears[length(allYears)])
yearsSeq <- firstYear:lastYear

# Combine the years, intensity and duration objects, for ease of comparison
data.frame(years=yearsSeq, 
           durationValue=locationDuration[,1], 
           intensityValue=locationIntensity[,1])

The output now closely resembles the output format that would be achieved if the location were extracted prior to analysis. Pulling the remaining stats (“firstMaxValue”, “secondMaxValue”, and “min”) and combining those as well would result in an identical result.