knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.width = 8, fig.height = 5 ) NOT_CRAN <- identical(tolower(Sys.getenv("NOT_CRAN")), "true") knitr::opts_chunk$set(purl = NOT_CRAN, eval = NOT_CRAN)
The eddi R package provides easy access to the Evaporative Demand Drought Index (EDDI) data - an experimental drought monitoring and early warning guidance tool produced by the National Oceanic and Atmospheric Administration. EDDI is available at multiple timescales, from weekly to monthly, providing insight into short-term flash droughts, and long-term droughts. More information on the EDDI product is available on the NOAA EDDI homepage.
This vignette covers a common use case for EDDI: acquiring data over a region of interest defined by a shapefile, masking the EDDI data to that region, and saving GeoTIFF files containing EDDI data for the region of interest.
By default, the eddi package returns data for the continental United States, southern parks of Canada, and northern parts of Mexico. But, you may only be interested in a region of interest, as defined by a shapefile. Here, you will load a shapefile for the state of North Carolina that is distributed by default with the sf package.
library(sf) library(raster) library(eddi) roi <- st_read(system.file("shape/nc.shp", package="sf"))
If you are using a different shapefile, replace
system.file("shape/nc.shp", package="sf") with its file path, e.g.,
roi object contains multiple columns of data, and a
that contains spatial information on the region of interest, which in this
case consists of multiple counties.
Because you don't necessarily care about each county, but rather you want the entire state (including all counties) you can use a spatial union to join data from all counties:
roi <- st_union(roi) roi
To acquire EDDI data, you can use the
You will fetch the 1 week timescale data for July, 2018:
eddi_raster <- get_eddi(date = "2018-07-01", timescale = "1 week")
eddi_raster object is a
RasterStack with one layer, and you can see
information on the spatial extent, resolution, and coordinate reference
system by printing the object:
Plot the data with a custom color palette to see what the data look like:
color_pal <- colorRampPalette(c("blue", "lightblue", "white", "pink", "red")) plot(eddi_raster, col = color_pal(255))
Now you want to subset or mask the EDDI data to the region of interest. First, you need to ensure that the raster data and the polygon for the region of interest have the same coordinate reference system.
roi_reprojected <- st_transform(roi, crs = projection(eddi_raster))
Now, graphically verify that they align as expected:
plot(eddi_raster, col = color_pal(255)) plot(roi_reprojected, add = TRUE)
Now, you can crop the EDDI data to match extents with the region of interest,
then mask the raster set all values outside of the region of interest to
Because the raster package requires sp objects, rather than sf objects, you
will coerce our roi to a sp object first.
roi_sp <- as(roi_reprojected, 'Spatial') cropped_eddi <- crop(eddi_raster, roi_sp) masked_eddi <- mask(cropped_eddi, roi_sp)
You can plot the masked EDDI raster along with the ROI to confirm:
plot(masked_eddi, col = color_pal(255)) plot(roi_sp, add = TRUE)
To write a GeoTIFF file of our
masked_eddi object, you can use
You can modify the
output_directory below to save this file in a particular
location on your filesystem.
output_directory <- tempdir() output_file <- file.path(output_directory, 'eddi-over-roi.tif') writeRaster(masked_eddi, output_file)
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.