add OpenET download function and comparison script

Author: ashiklom

modules/data.land/NAMESPACE

@@ -29,6 +29,7 @@ export(extract_SM_CDS)
 export(extract_chirps_remote)
 export(extract_cimis_date)
 export(extract_cimis_dates)
+export(extract_openet_daily)
 export(extract_soil_gssurgo)
 export(extract_soil_nc)
 export(extract_veg)

modules/data.land/R/openET.R

@@ -0,0 +1,82 @@
+#' Extract daily ET data from OpenET
+#'
+#' Note that this requires the environment variable `OPENET_API_KEY` to be set.
+#' A convenient way to do this is via a `.Renviron`, either globally
+#' (`~/.Renviron`) or in the current working directory (`./.Renviron`), with
+#' contents like:
+#'
+#' ```
+#' OPENET_API_KEY="abcdefg123456"
+#' ```
+#'
+#' You can obtain an OpenET API key from the OpenET data portal.
+#'
+#' @param design_points `data.frame` of design points with columns `lat` and `lon`
+#' @param start_date Start date for data extraction
+#' @param end_date End date for data extraction
+#'
+#' @return `design_points` `data.frame` with additional columns `date`, and
+#' `et_mm_day` (ET, mm/day)
+#' @export
+extract_openet_daily <- function(design_points, start_date, end_date) {
+  api_key <- Sys.getenv("OPENET_API_KEY")
+  if (api_key == "") {
+    stop("OPENET_API_KEY environment variable is not set")
+  }
+
+  start_date_str <- format(start_date, "%Y-%m-%d")
+  end_date_str <- format(end_date, "%Y-%m-%d")
+
+  request_body_template <- list(
+    date_range = c(start_date_str, end_date_str),
+    interval = "daily",
+    model = "Ensemble",
+    variable = "ET",
+    reference_et = "gridMET",
+    units = "mm",
+    file_format = "JSON"
+  )
+
+  prep_request <- function(lon, lat) {
+    request_body <- request_body_template
+    request_body$geometry <- c(lon, lat)
+
+    httr2::request("https://openet-api.org/raster/timeseries/point") |>
+      httr2::req_headers(Authorization = api_key) |>
+      httr2::req_body_json(request_body) |>
+      httr2::req_throttle(capacity = 10, fill_time_s = 1) |>
+      httr2::req_retry(max_tries = 3) |>
+      httr2::req_timeout(seconds = 150)
+  }
+
+  raw_results <- design_points |>
+    dplyr::mutate(
+      reqs = purrr::map2(.data$lon, .data$lat, prep_request),
+      resps = httr2::req_perform_parallel(
+        raw_results[["reqs"]],
+        max_active = 10,
+        on_error = "continue"
+      )
+    )
+
+  parse_response <- function(resp) {
+    if (!inherits(resp, "httr2_response")) {
+      return(NULL)
+    }
+    data <- httr2::resp_body_json(resp)
+    if (length(data) == 0 || is.null(data[[1]]$time)) {
+      return(NULL)
+    }
+    tibble::tibble(
+      date = as.Date(purrr::map_chr(data, "time")),
+      et_mm_day = purrr::map_dbl(data, "et")
+    )
+  }
+
+  results <- raw_results |>
+    dplyr::mutate(results = purrr::map(.data$resps, parse_response)) |>
+    dplyr::select(-c("reqs", "resps")) |>
+    tidyr::unnest("results", keep_empty = TRUE)
+
+  results
+}

modules/data.land/inst/compare-cimis-openet.R

@@ -0,0 +1,41 @@
+#!/usr/bin/env Rscript
+
+if (interactive()) {
+  devtools::load_all("modules/data.land")
+} else {
+  library(PEcAn.data.land)
+}
+
+design_points_path <- "~/projects/cimis-to-irrigation/design_points.csv"
+cimis_eto_cog_path <- "~/data/CIMIS-ETo-COG"
+
+dates <- seq.Date(as.Date("2020-03-01"), as.Date("2020-11-30"), "day")
+design_points <- readr::read_csv(design_points_path) |>
+  head(10)
+
+cimis_et <- extract_cimis_dates(
+  design_points,
+  dates,
+  cimis_eto_cog_path,
+  .progress = TRUE
+)
+
+openet_et <- extract_openet_daily(
+  design_points,
+  min(dates),
+  max(dates)
+)
+
+combined <- cimis_et |>
+  dplyr::rename(cimis = etref_mm_day) |>
+  dplyr::full_join(openet_et |> dplyr::rename(openet = et_mm_day))
+
+combined_long <- combined |>
+  tidyr::pivot_longer(c(cimis, openet), names_to = "source", values_to = "et_mm_day")
+
+library(ggplot2)
+ggplot(combined_long) +
+  aes(x = date, y = et_mm_day, color = source) +
+  geom_line() +
+  facet_wrap(~id) +
+  theme_bw()