CONUS Drought Monitoring Pipeline: Methods

Zachary H. Hoylman, R. Kyle Bocinsky, and Kelsey G. Jencso — Montana Climate Office, University of Montana

Updated: Version 1.2

Abstract The Montana Climate Office (MCO) CONUS Drought Monitoring Pipeline is an operational system that computes a suite of standardized drought indices on a daily basis across the contiguous United States (CONUS). The pipeline ingests daily gridded meteorological fields from gridMET [1], a high-resolution (~4 km; 1/24°) observationally constrained climate dataset that merges PRISM spatial patterns with reanalysis-derived temporal signals. Input variables include precipitation, potential evapotranspiration (PET), vapor pressure deficit (VPD), and maximum air temperature.

From these inputs the pipeline derives four standardized drought index families — the Standardized Precipitation Index (SPI), the Standardized Precipitation-Evapotranspiration Index (SPEI), the Evaporative Demand Drought Index (EDDI), and the Standardized VPD Index (SVPDI) — together with maximum temperature percentiles (Tmax) and ancillary products including precipitation accumulation (mm), percent of normal and deviation from normal for both precipitation and VPD, empirical percentile ranks for precipitation, VPD, and Tmax, and temperature deviation from normal. Each metric is computed across multiple aggregation timescales ranging from 15 to 730 days, as well as water-year and year-to-date accumulations.

A distinguishing feature of this pipeline is its support for flexible climatological reference periods. Because a warming climate produces non-stationary baseline statistics, drought classifications derived from a fixed historical normal can be misleading about contemporary conditions [2]. The pipeline therefore supports rolling 30-year windows, fixed WMO-standard normals, and full period-of-record baselines, all configurable at runtime. All outputs are written as Cloud-Optimized GeoTIFFs (COG) with LZW compression, tiled internally for efficient partial reads by downstream visualization and analysis tools.

2. Data

2.1 gridMET Climate Data

All meteorological inputs are obtained from gridMET [1], a daily surface meteorological dataset covering CONUS at a spatial resolution of approximately 4 km (1/24°). gridMET blends the temporal attributes of regional-scale reanalysis fields (NLDAS-2) with the spatial climatological patterns embedded in PRISM (Parameter-elevation Regressions on Independent Slopes Model), yielding a dataset that is both temporally consistent and spatially detailed. The archive extends from 1979 to near-present (updated daily), providing the multi-decadal record required for robust climatological fitting. Data are freely available at climatologylab.org/gridmet.html.

The following gridMET variables are used by this pipeline:

Variable gridMET name Units Used for
Precipitationprmm/daySPI, SPEI, precip % of normal, precip deviation, precip percentile, precip accumulation (mm)
Reference ET (ETo)petmm/daySPEI, EDDI
Vapor pressure deficitvpdkPaSVPDI, VPD % of normal, VPD deviation, VPD percentile
Maximum air temperaturetmmxKTmax percentile, Tmax deviation

Download & Caching Strategy

gridMET distributes data as annual NetCDF-4 files. The pipeline loops over every year from START_YEAR (default 1979) through the current calendar year for all four variables (pr, pet, vpd, tmmx), using timestamp-conditional downloads via curl -R -z:

  • -R — Sets the local file's modification time to the server's Last-Modified timestamp on download.
  • -z <file> — Sends an If-Modified-Since header. If the remote file has not been updated since the local copy's mtime, the download is skipped entirely (HTTP 304).

This means every year — including historical years — is checked on each run, but only files that have actually changed on the server are re-downloaded. In practice, completed historical years are rarely re-fetched after the initial download (only when the provider updates or corrects historical fields), while the current (in-progress) year file is updated whenever the provider appends new daily observations. Failed downloads are retried up to 3 times with a short backoff.

All annual NetCDFs are retained to allow reprocessing with alternative reference periods without re-downloading. The default archive start year is 1979 (the full gridMET record); this is configurable via the START_YEAR environment variable.

3. Climatological Reference Periods

Scientific Motivation

Standardized drought indices measure anomalies relative to a baseline climatology. The choice of baseline profoundly affects index values: for example, the same temperature, reference ET, or precipitation value could indicate moderate drought against a 1950–1980 baseline but near-normal conditions against a 1991–2020 WMO normal. In other words, as the climate warms, the statistical properties of temperature, PET, and water balance may trend systematically in any of the distribution’s moments (mean, standard deviation, skewness, kurtosis, etc.), violating the stationarity assumption embedded in fixed-period normals. Hoylman et al. (2022) [2] demonstrate that non-stationarity in baseline climate substantially alters drought frequency and severity assessments, motivating the use of flexible, contemporary reference periods for operational monitoring.

Supported Modes

The reference period is controlled at runtime via the CLIM_PERIODS environment variable (comma-separated, multiple periods may be computed in a single run):

Spec syntax Mode Slug appended to filename Description
rolling:30 Rolling 30-year window rolling-30 Rolling 30-year window ending at the current year; tracks contemporary climatology (default). rolling:N generalizes to any window length.
fixed:1991:2020 Fixed interval fixed-1991-2020 WMO 1991–2020 standard normal; suitable for comparison with published climatologies
full Full period of record full All years from START_YEAR to present; maximizes sample size

Example configurations:

# Single rolling 30-year window (default)
CLIM_PERIODS="rolling:30"

# Two periods computed in one run
CLIM_PERIODS="rolling:30,fixed:1991:2020"

# Fixed period only
CLIM_PERIODS="fixed:1991:2020"

# Full period of record
CLIM_PERIODS="full"

Each period produces a fully independent set of output files with the slug appended before the .tif extension. Intermediate tile directories are also slug-tagged to prevent collisions when multiple periods are computed simultaneously.

Minimum anchor requirement: At least 10 calendar-year anchors must fall within the selected reference period for any distribution fit to proceed at a given grid cell. Cells that do not meet this threshold are assigned NA for that index and timestep.

4. Drought Index Methods

All indices are computed at each grid cell independently. For a given aggregation window of W days, a rolling sum or mean of the input variable is first computed over the full historical record. The resulting time series is then stratified by calendar day-of-year, and a probability distribution (or empirical CDF) is fitted to the subsample of values that fall within the selected climatological reference period. The current (most recent) observation is then transformed to a standardized anomaly or probability using that fitted distribution.

Note on aggregation: For fixed-length windows, the choice of rolling sum vs. rolling mean is invariant under the subsequent standardization — the two differ by a constant scale factor (1/W) that cancels in the fitted CDF, so the resulting standardized indices are identical regardless of which is used. The operational pipeline uses rolling sums for SPI, SPEI, and EDDI (yielding outputs in physically interpretable accumulated units) and rolling means for SVPDI (so VPD is reported in its native units).

Sign conventions. Two distinct sign conventions are used across the indices on this page, reflecting the underlying variable. Hydroclimatic standardized indices (SPI, SPEI) follow the precipitation-anomaly convention: negative values indicate drought, positive values indicate wet anomalies. Atmospheric-demand indices (EDDI, SVPDI) follow the demand-anomaly convention: positive values indicate drier-than-normal atmospheric conditions, negative values indicate below-normal demand. The Tmax percentile is not a centered z-score but a probability metric where higher values indicate unusually warm conditions.

4.1 Standardized Precipitation Index (SPI)

Motivation

The SPI [3] expresses precipitation anomalies in units of standard deviations relative to a fitted climatological distribution, making it comparable across locations and timescales with differing precipitation regimes. It is the most widely used standardized drought index and serves as a fundamental building block for multi-variable composite assessments.

Input

Daily precipitation (pr), in millimeters.

Aggregation

Rolling window sums over the following timescales:

  • Fixed windows: 15, 30, 45, 60, 90, 120, 180, 365, and 730 days
  • Water-year accumulation (Oct 1 – current date, or equivalent)
  • Year-to-date accumulation (Jan 1 – current date)

Statistical Method

Precipitation is inherently bounded at zero and right-skewed; a two-parameter Gamma distribution (shape α, scale β) is therefore fitted to the climatological sample for each grid cell and aggregation window. Parameter estimation follows the L-moment / probability-weighted moment (PWM) approach of Hosking (1990) [4], as implemented in the R package lmomco:

  1. Zero handling (Stagge et al., 2015) — Following the methodology of Stagge et al. (2015) [5], zero (or below-threshold) precipitation values are handled via a mixed distribution approach using the Weibull plotting position. The probability of zero precipitation is estimated as p0 = m / (n + 1), where m is the number of zero values and n is the total sample size. Zero-valued observations are assigned the center of probability mass below the threshold:
    p̄0 = (m + 1) / [ 2(n + 1) ]
    Non-zero values are mapped to the upper portion of the probability space:
    H(x) = p0 + (1 − p0) ⋅ G(x; α, β)
    where G is the Gamma CDF fitted only to positive values. This approach correctly represents the discrete mass at zero without the artifacts introduced by replacing zeros with small positive constants.
  2. PWM estimation — Unbiased probability-weighted moments β0, β1 are computed from the positive (non-zero) values in the climatological sample.
  3. L-moment conversion — PWMs are converted to L-moments (λ1, λ2, τ3).
  4. Parameter estimation — lmomco::pargam() returns (α, β) from L-moments of the positive values only.
  5. CDF evaluation — For non-zero observations, the mixed CDF H(x) is evaluated to obtain a non-exceedance probability p. For zero-valued observations, p = p̄0 (center of probability mass).
  6. Normal quantile transform — SPI = Φ−1(p), where Φ−1 is the standard normal quantile function. No clamping is applied to the resulting SPI values; in the rare case that the fitted CDF saturates numerically (p = 0 or 1), the resulting non-finite value is set to NA.
SPI = Φ−1[ H( XW ) ]   (exact zeros take p = p̄0)

Interpretation

Category Description Percentile Range SPI / SPEI Values
Normal or wet conditionsx > 30x > −0.524
D0Abnormally Dry30 ≥ x > 20−0.524 ≥ x > −0.842
D1Moderate Drought20 ≥ x > 10−0.842 ≥ x > −1.282
D2Severe Drought10 ≥ x > 5−1.282 ≥ x > −1.645
D3Extreme Drought5 ≥ x > 2−1.645 ≥ x > −2.054
D4Exceptional Droughtx < 2x < −2.054

Minimum Data Requirement

At least 10 anchor years must fall within the selected reference period for a distribution fit to be attempted. Grid cells with fewer anchor years receive a fill value of NA for the current timestep. In addition, the Gamma fit itself requires at least 3 positive (non-zero) values with non-zero variance in the reference sample; cells that fail this condition (e.g. a reference period that is nearly all zeros at short timescales in arid regions) also return NA.

4.2 Ancillary Precipitation Metrics

In addition to the primary drought indices, four ancillary metrics are derived directly from rolling precipitation sums to provide operational context:

Metric Formula Units Description
Percent of Normal (PON) ( XW / μclim ) × 100 % Current accumulation as a percentage of the climatological mean
Departure XW − μclim mm Absolute anomaly relative to climatological mean
Percentile F̂n( XW ) 0–1 Empirical ECDF rank within the climatological sample
Accumulation (precip-mm) XW mm Raw rolling precipitation sum

The climatological mean (μclim) and ECDF are derived from the same reference period used for SPI fitting, ensuring internal consistency across all precipitation-based products.

4.3 Standardized Precipitation-Evapotranspiration Index (SPEI)

Motivation

The SPEI [6] extends SPI by incorporating the atmospheric evaporative demand, making it sensitive to warming-driven increases in PET even when precipitation is unchanged. In a warming climate, SPEI captures atmospheric-demand-driven drought stress that SPI would miss.

Input

Water balance = precipitation (pr) − reference evapotranspiration (pet), in mm/day. Negative values indicate a moisture deficit; positive values indicate a moisture surplus.

Aggregation

Rolling window sums of the daily water balance over the same timescales as SPI; equivalently, the rolling sum of precipitation minus the rolling sum of PET, both accumulated over the same window. Units are millimeters of accumulated moisture surplus (positive) or deficit (negative). The Generalized Logistic (GLO) distribution is then fitted to this aggregated quantity for each calendar day-of-year within the climatological reference period, and the current observation is transformed via the same fitted distribution.

Statistical Method

The water balance variable is unbounded below (large PET on dry days) and has heavier distributional tails than precipitation alone. The GLO distribution, which accommodates both heavy tails and asymmetry, is therefore used in place of the Gamma:

  1. L-moments estimated from the climatological sample via unbiased PWMs.
  2. lmomco::parglo() fits the three-parameter GLO (ξ, α, k).
  3. GLO CDF evaluated at current observation → probability p.
  4. SPEI = Φ−1(p).
SPEI = Φ−1[ FGLO(ξ,α,k)( P − PET )W ]

Interpretation

SPEI values follow the same classification thresholds as SPI. Because it incorporates evaporative demand, SPEI typically shows stronger drying trends over recent decades than SPI alone, consistent with the non-stationarity findings of [2].

4.4 Evaporative Demand Drought Index (EDDI)

Motivation

EDDI [7] isolates the atmospheric demand component of the hydrological cycle. Because it is derived solely from PET — without precipitation — it captures early warning signals of drought stress driven by high temperature, low humidity, or high wind speed before soil moisture depletion has occurred.

Input

Reference evapotranspiration (pet), in mm/day.

Aggregation

Rolling window sums (cumulative evaporative demand) over the same timescales as SPI.

Statistical Method

EDDI uses a nonparametric rank-based approach, avoiding assumptions about the distributional form of PET:

  1. Rank — Observations within the climatological sample are ranked from highest (rank 1) to lowest. A rank of 1 therefore corresponds to maximum evaporative demand (drought stress).
  2. Tukey plotting position — Following Hobbins et al. (2016) [7], the empirical probability is estimated via the Tukey plotting position:
    pi = (i − 0.33) / (n + 0.33)
    where i is the rank of the aggregated E0 in the historical time series (i = 1 for maximum E0) and n is the number of observations being ranked. When the reference period is rolling or full-record, the current observation is the final member of the climatological sample (n = sample size); under a fixed reference period that excludes the current year, the current observation is appended to the sample and ranked within the resulting n + 1 values. Because rank 1 = highest demand, p is small for drought conditions.
  3. Inverse normal approximation — EDDI is derived via the rational approximation of Abramowitz and Stegun (1965) [8], following the formulation in Hobbins et al. (2016) [7] and Vicente-Serrano et al. (2010) [10]. For p ≤ 0.5, W = √(−2⋅ln(p)); for p > 0.5, p is replaced with 1−p and the sign of EDDI is reversed:
    EDDI = W − ( C0 + C1W + C2W2 ) / ( 1 + d1W + d2W2 + d3W3 )
    with constants C0 = 2.515517, C1 = 0.802853, C2 = 0.010328, d1 = 1.432788, d2 = 0.189269, and d3 = 0.001308, as tabulated in Hobbins et al. (2016).
  4. Sign convention — High PET (low p) produces positive EDDI values, indicating drier-than-normal conditions. Negative values indicate wet anomalies. A zero EDDI indicates that accumulated E0 equals the climatological median.
EDDI = −Φ−1[ pTukey( E0,W ) ]   (rank 1 = highest E0; the leading minus enforces positive EDDI for drought)

Interpretation

Positive EDDI indicates above-normal atmospheric demand (drought precursor or ongoing drought amplifier). EDDI provides a pure atmospheric signal unconfounded by antecedent precipitation and is particularly valuable for flash drought early warning.

4.5 Standardized VPD Index (SVPDI)

Motivation

Vapor pressure deficit (VPD) is the difference between the saturation vapor pressure and the actual vapor pressure. Elevated VPD increases transpiration demand and can cause stomatal closure even when soil moisture is adequate, directly linking atmospheric drought to vegetation stress. The SVPDI [9] provides a standardized measure of VPD anomaly that is complementary to the precipitation-based indices.

Input

Mean daily vapor pressure deficit (vpd), in kPa.

Aggregation

Rolling window means over the same timescales as SPI.

Statistical Method

Following Nwayor et al. (2024) [9], VPD is treated as a zero-limited, right-skewed variable and fitted with a Gamma distribution. Nwayor et al. compared Gamma, Weibull, and log-normal distributions via AIC and found the Gamma to be the best-performing distribution for VPD across most locations globally. Our implementation uses L-moment/PWM parameter estimation (rather than MLE) for consistency with the SPI fitting procedure, and applies the same Stagge et al. (2015) [5] mixed-distribution zero handling used for SPI:

  1. Zero handling (Stagge et al., 2015) — Zero VPD values (rare but possible at high latitudes or during cold conditions) are handled via the same mixed-distribution approach as SPI. The probability of zero is estimated as p0 = m / (n + 1), and zero-valued observations are assigned the center of probability mass: p̄0 = (m + 1) / [2(n + 1)]. The Gamma distribution is fitted only to positive values.
  2. Gamma fitting — Unbiased PWMs → L-moments → lmomco::pargam() on positive values only.
  3. CDF evaluation — For non-zero observations: H(x) = p0 + (1 − p0) ⋅ G(x; α, β). For zero observations: p = p̄0.
  4. Normal quantile transform — SVPDI = Φ−1(p). No clamping is applied.
SVPDI = Φ−1[ H( VPDW ) ]

Interpretation

Positive SVPDI indicates above-normal VPD (heightened atmospheric demand and potential vegetation stress). Negative SVPDI indicates below-normal VPD (reduced atmospheric demand). SVPDI follows the atmospheric-demand sign convention shared with EDDI — positive values indicate drier-than-normal atmospheric conditions — which is the opposite of the SPI/SPEI hydroclimatic convention, where negative values indicate drought.

4.6 Ancillary VPD Metrics

Alongside the SVPDI, three ancillary metrics are computed from the rolling VPD means to provide operational context:

Metric Formula Units Description
Percent of Normal (PON) ( VPDW / μclim ) × 100 % Current VPD as a percentage of the climatological mean
Departure VPDW − μclim kPa Absolute anomaly relative to climatological mean VPD
Percentile F̂n( VPDW ) 0–1 Empirical ECDF rank within the climatological sample

4.7 Maximum Temperature Metrics (Tmax)

Motivation

Anomalously high temperatures amplify drought through increased PET and accelerated snowmelt. A percentile-based temperature index provides context for heat as a drought driver without coupling it to a precipitation signal.

Input

Daily maximum air temperature (tmmx), in Kelvin.

Aggregation

Rolling window means over the standard timescales.

Statistical Method

A nonparametric empirical CDF approach is used; no parametric distribution is assumed:

  1. The empirical cumulative distribution function (ECDF) is constructed from the climatological sample: F̂n(x) = (number of reference values ≤ x) / n (evaluated via a vectorized rank computation equivalent to base R ecdf()).
  2. The ECDF is evaluated at the current observation to produce a probability (0–1 scale), where 1 indicates the warmest value on record within the reference period.
Tmaxpctile = F̂n( TmaxW )

Deviation from Normal

In addition to the percentile, the pipeline computes the deviation of the current rolling-mean Tmax from the climatological mean for each timescale:

Tmaxdev = TmaxW − μclim

Positive values indicate above-normal temperatures; negative values indicate below-normal. Units are K.

Interpretation

For the Tmax percentile: values > 0.9 indicate exceptionally warm conditions (top decile of the reference climatology). The percentile is reported as a 0–1 probability rather than a standardized anomaly, making it directly interpretable as an exceedance percentile. For the Tmax deviation: values are the absolute departure of the current rolling-mean maximum temperature from the climatological mean, in K (equivalent to °C for a temperature difference); positive values indicate above-normal temperatures.

5. Spatial Processing Architecture

Tiling

Loading the full 47-year daily record for all of CONUS at once would require far more memory than is practical, so the domain is partitioned into rectangular 1°×1° tiles (configurable via the TILE_DX/TILE_DY environment variables) and processed tile-by-tile. This keeps the per-worker memory footprint small enough to run on modestly sized (cost-effective) ECS Fargate tasks while still exploiting all available cores. Each tile spans approximately 110×85 km at mid-latitudes, containing roughly 600 grid cells (24×24). Tiles are processed independently, which enables embarrassingly parallel execution and constrains peak memory usage per process. Tiles that fail (e.g. due to memory pressure) are automatically subdivided into four quadrants and retried at the smaller size.

Parallelism

Within each R script, tiles are distributed across CPU cores using pbmcapply::pbmclapply(), a fork-based parallel backend (equivalent to parallel::mclapply()) that additionally renders a real-time text progress bar. Key design constraints:

Mosaicking

After all tiles for a given index and timescale are written, the pipeline assembles them into a seamless CONUS mosaic:

  1. A GDAL Virtual Raster (VRT) is constructed referencing all tile GeoTIFF files. VRT creation is instantaneous (no pixel copying) and is used to drive the subsequent translation.
  2. gdal_translate converts the VRT to a Cloud-Optimized GeoTIFF (-of COG) with LZW compression (COMPRESS=LZW, PREDICTOR=2). The COG driver applies internal 512×512 tiling and automatically generates the overview pyramid, ensuring efficient partial reads and thumbnail rendering by web map clients.

Land Mask

An optional CONUS land mask derived from Natural Earth state boundary polygons can be applied to suppress ocean and non-CONUS pixels in the final mosaic. This is controlled by the CONUS_MASK environment variable.

6. Output Products

All final output layers are written to the conus_drought/ directory as Cloud-Optimized GeoTIFFs. Filenames follow the pattern: 

{index}_{timescale}_{clim_slug}_{YYYY-MM-DD}.tif
# e.g.: spi_30d_rolling-30_2026-03-10.tif
#       spei_90d_fixed-1991-2020_2026-03-10.tif
#       eddi_365d_full_2026-03-10.tif

where the trailing date is the data date (the last day of the aggregation window). Published copies in the latest/ directory carry the same names with the date suffix stripped (e.g. spi_30d_rolling-30.tif), so operational consumers can request a stable URL; the accompanying manifest.csv and latest-date.txt record the data date of each layer.

The {clim_slug} component identifies the climatological reference period the layer was standardized against (Section 3): rolling-30 = rolling 30-year window ending at the current year; fixed-1991-2020 = the fixed WMO 1991–2020 standard normal; full = the full period of record (START_YEAR to present).

The units listed below describe the gridded data as written by this pipeline — i.e. what is stored in the public S3 bucket at data2.climate.umt.edu/gridmet/derived/conus_drought/. Downstream applications such as the D³ dashboard may convert units on the fly for display (e.g. mm to inches), but the underlying rasters always carry the units below.

Index Timescales Units / Range Distribution
SPI 15d, 30d, 45d, 60d, 90d, 120d, 180d, 365d, 730d, WY, YTD Dimensionless z-score Gamma (L-moments)
SPEI 15d, 30d, 45d, 60d, 90d, 120d, 180d, 365d, 730d, WY, YTD Dimensionless z-score Generalized Logistic (L-moments)
EDDI 15d, 30d, 45d, 60d, 90d, 120d, 180d, 365d, 730d, WY, YTD Dimensionless z-score Nonparametric (Tukey plotting position)
SVPDI 15d, 30d, 45d, 60d, 90d, 120d, 180d, 365d, 730d, WY, YTD Dimensionless z-score Gamma (L-moments)
Tmax percentile 15d, 30d, 45d, 60d, 90d, 120d, 180d, 365d, 730d, WY, YTD Probability [0, 1] Empirical ECDF
Percent of Normal Same as SPI % (unbounded) Climatological mean
Departure Same as SPI mm Climatological mean
Percentile (precip) Same as SPI Probability [0, 1] Empirical ECDF
Accumulation (precip-mm) Same as SPI mm Raw sum
VPD % of Normal Same as SVPDI % (unbounded) Climatological mean
VPD Departure Same as SVPDI kPa Climatological mean
VPD Percentile Same as SVPDI Probability [0, 1] Empirical ECDF
Tmax deviation Same as Tmax percentile K Climatological mean

All rasters share the native gridMET CRS (EPSG:4326, geographic/WGS84) and pixel size (1/24° ≈ 4 km). In the primary (analysis-grade) products, values are stored as 32-bit floats and NoData is encoded as NaN (IEEE 754 32-bit float).

Web-Optimized Copies

Alongside the primary products, the pipeline publishes a parallel set of web-optimized COGs (conus_drought_web/) intended for browser-based visualization. Values are multiplied by 100 and stored as integers (Int16; Int32 for wide-range metrics such as precipitation accumulation and percent of normal), with SCALE=0.01 recorded in the GeoTIFF metadata so GDAL-aware clients recover physical units automatically. These files use DEFLATE compression, include a full overview pyramid, and encode NoData as −9999. The primary Float32 products remain the authoritative source for analysis.

7. Software & Reproducibility

Runtime Environment

The pipeline runs inside a Docker container based on the rocker/geospatial image (R 4.4.1, GDAL 3.x, PROJ 9.x, GEOS). All system libraries are pinned at image-build time, ensuring byte-for-byte reproducibility of spatial operations.

Key R Packages

PackageRole
terraRaster I/O, spatial operations, masking
sfVector geometry (tile boundaries, land mask)
lmomcoL-moment parameter estimation; Gamma, GLO distributions
ncdf4Reading gridMET NetCDF-4 source files
fsFile system operations (path manipulation, directory management)
purrrFunctional programming utilities for list operations
pbmcapplyFork-based parallel apply with real-time progress bar
gdalUtilitiesR wrappers for gdal_translate, gdalbuildvrt

Configurable Parameters (Environment Variables)

VariableDefaultDescription
START_YEAR1979First year of gridMET archive to download/use
CLIM_PERIODSrolling:30Comma-separated reference period specification(s)
CORES12Number of parallel worker processes
TILE_DX1Tile width in degrees
TILE_DY1Tile height in degrees
CONUS_MASK1Apply CONUS land mask (0 = no mask)
DATA_DIR~/mco-drought-conus-dataRoot directory for input NetCDF cache and outputs
TIMESCALES15,30,45,60,90,120,180,365,730,wy,ytdTimescales to compute
TILE_IDS(unset)Subset of tile IDs to process; unset = all tiles

Docker Usage

# Build the image
docker compose build

# Run the full pipeline with default settings
docker compose up

# Run with multiple reference periods
docker compose run --rm -e CLIM_PERIODS="rolling:30,fixed:1991:2020" mco-drought

# Run with limited cores for debugging
docker compose run --rm -e CORES=4 mco-drought

# Quick test (single metric, few tiles)
docker compose run --rm mco-drought bash pipeline/run_test.sh

Source Code

Source code and Docker configuration are maintained in the mco-drought-conus repository. The pipeline is structured as follows: 

mco-drought-conus/
├── Dockerfile
├── docker-compose.yml
├── R/
│   ├── drought-functions.R        # Core statistical fitting functions
│   ├── pipeline-common.R          # Shared pipeline infrastructure
│   ├── 1_gridmet-cache.R          # gridMET download and caching
│   ├── 2_metrics-precip.R         # SPI + ancillary metrics
│   ├── 3_metrics-spei.R           # SPEI
│   ├── 4_metrics-eddi.R           # EDDI
│   ├── 5_metrics-vpd.R            # SVPDI
│   └── 6_metrics-tmax.R           # Tmax percentile + deviation
├── pipeline/
│   ├── run_once.sh                # Orchestration shell script (container entry point)
│   ├── make_web_cogs.sh           # Web-optimized COG post-processing
│   └── run_test.sh                # Quick-test runner
├── scripts/                       # Deployment utilities (ECR push, etc.)
├── terraform/                     # AWS infrastructure as code
└── docs/
    └── methods.html               # This document

8. Partners & Acknowledgments

The Drought Data Dashboard (D³) is a collaboration between the following institutions. This work is supported by ongoing partnerships across academic, state, and federal drought monitoring programs.

Upstream data providers and authoritative sources (gridMET, NOHRSC SNODAS, NOAA GHCNd, USGS NWIS, U.S. Drought Monitor, BIA, USGS WBD) are cited in References and on each dataset's inline tooltip within the dashboard.

9. References

  1. Abatzoglou JT. Development of gridded surface meteorological data for ecological applications and modelling. International Journal of Climatology. 2013;33(1):121–131. doi:10.1002/joc.3413
  2. Hoylman ZH, Bocinsky RK, Jencso KG. Drought assessment has been outpaced by climate change: empirical arguments for a paradigm shift. Nature Communications. 2022;13:2715. doi:10.1038/s41467-022-30316-5
  3. McKee TB, Doesken NJ, Kleist J. The relationship of drought frequency and duration to time scales. Proceedings of the 8th Conference on Applied Climatology. 1993;17(22):179–183. American Meteorological Society.
  4. Hosking JRM. L-moments: analysis and estimation of distributions using linear combinations of order statistics. Journal of the Royal Statistical Society, Series B. 1990;52(1):105–124. doi:10.1111/j.2517-6161.1990.tb01775.x
  5. Stagge JH, Tallaksen LM, Gudmundsson L, Van Loon AF, Stahl K. Candidate distributions for climatological drought indices (SPI and SPEI). International Journal of Climatology. 2015;35(13):4027–4040. doi:10.1002/joc.4267
  6. Beguería S, Vicente-Serrano SM, Reig F, Latorre B. Standardized precipitation evapotranspiration index (SPEI) revisited: parameter fitting, evapotranspiration models, tools, datasets and related applications. International Journal of Climatology. 2014;34(10):3001–3023. doi:10.1002/joc.3887
  7. Hobbins MT, Wood A, McEvoy DJ, Huntington JL, Morton C, Anderson M, Hain C. The Evaporative Demand Drought Index. Part I: Linking drought evolution to variations in evaporative demand. Journal of Hydrometeorology. 2016;17(6):1745–1761. doi:10.1175/JHM-D-15-0121.1
  8. Abramowitz M, Stegun IA, eds. Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables. National Bureau of Standards Applied Mathematics Series 55. Washington, DC: U.S. Government Printing Office; 1965.
  9. Nwayor IJ, Robeson SM, Ficklin DL, Maxwell JT. A multiscalar standardized vapor pressure deficit index for drought monitoring and impacts. International Journal of Climatology. 2024;44:5825–5838. doi:10.1002/joc.8668
  10. Vicente-Serrano SM, Beguería S, López-Moreno JI. A multiscalar drought index sensitive to global warming: the Standardized Precipitation Evapotranspiration Index. Journal of Climate. 2010;23(7):1696–1718. doi:10.1175/2009JCLI2909.1