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.,
st_read("path/to/file.shp")
.
The roi
object contains multiple columns of data, and a geometry
column
that contains spatial information on the region of interest, which in this
case consists of multiple counties.
roi
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 get_eddi()
function.
You will fetch the 1 week timescale data for July, 2018:
eddi_raster <- get_eddi(date = "2018-07-01", timescale = "1 week")
The 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:
eddi_raster
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 NA
.
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 writeRaster
.
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.