Get EDDI data over a region of interest and save a GeoTIFF

Maxwell B. Joseph

2019-05-23

Source: vignettes/eddi-roi-tutorial.Rmd

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.

Defining a 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.

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.

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:

Acquiring EDDI data

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:

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))

Masking to the region of interest

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.

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)

Saving GeoTIFF output

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)