From e2d7e9b5fb63ea6e6f6c6768c916f39e1b8a2bca Mon Sep 17 00:00:00 2001
From: Mike Mahoney <mike.mahoney.218@gmail.com>
Date: Mon, 26 Feb 2024 10:48:26 -0500
Subject: [PATCH] Add get_naip_data()

Fixes #43
---
 NAMESPACE                           |  1 +
 NEWS.md                             |  4 +++
 R/get_stac_data.R                   | 42 +++++++++++++++++++++++++++++
 man/get_stac_data.Rd                | 25 +++++++++++++++++
 tests/testthat/test-get_stac_data.R | 18 +++++++++++++
 5 files changed, 90 insertions(+)

diff --git a/NAMESPACE b/NAMESPACE
index df3d438..b48365b 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -9,6 +9,7 @@ export(filter_bands)
 export(filter_platforms)
 export(get_dem)
 export(get_landsat_imagery)
+export(get_naip_imagery)
 export(get_sentinel1_imagery)
 export(get_sentinel2_imagery)
 export(get_stac_data)
diff --git a/NEWS.md b/NEWS.md
index 394a59a..f7a9c83 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,5 +1,9 @@
 # rsi (development version)
 
+* Added `get_naip_data()`, a function for getting National Agricultural Imagery
+  Program data from (by default) Planetary Computer. Data covers the continental
+  United States.
+
 * `stack_rasters()` will only rename bands if `band_names` is the same length as 
   the number of bands in the output raster (or missing, or defined by a 
   function). It will now warn you if these lengths are different. Previously, if 
diff --git a/R/get_stac_data.R b/R/get_stac_data.R
index 0d63717..5dee510 100644
--- a/R/get_stac_data.R
+++ b/R/get_stac_data.R
@@ -565,6 +565,48 @@ get_landsat_imagery <- function(aoi,
   do.call(get_stac_data, args)
 }
 
+#' @rdname get_stac_data
+#' @export
+get_naip_imagery <- function(aoi,
+                             start_date,
+                             end_date,
+                             ...,
+                             pixel_x_size = 1,
+                             pixel_y_size = 1,
+                             asset_names = "image",
+                             stac_source = "https://planetarycomputer.microsoft.com/api/stac/v1",
+                             collection = "naip",
+                             query_function = default_query_function,
+                             sign_function = sign_planetary_computer,
+                             rescale_bands = FALSE,
+                             output_filename = paste0(proceduralnames::make_english_names(1), ".tif"),
+                             composite_function = "merge",
+                             limit = 999,
+                             gdalwarp_options = c(
+                               "-r", "bilinear",
+                               "-multi",
+                               "-overwrite",
+                               "-co", "COMPRESS=DEFLATE",
+                               "-co", "PREDICTOR=2",
+                               "-co", "NUM_THREADS=ALL_CPUS"
+                             ),
+                             gdal_config_options = c(
+                               VSI_CACHE = "TRUE",
+                               GDAL_CACHEMAX = "30%",
+                               VSI_CACHE_SIZE = "10000000",
+                               GDAL_HTTP_MULTIPLEX = "YES",
+                               GDAL_INGESTED_BYTES_AT_OPEN = "32000",
+                               GDAL_DISABLE_READDIR_ON_OPEN = "EMPTY_DIR",
+                               GDAL_HTTP_VERSION = "2",
+                               GDAL_HTTP_MERGE_CONSECUTIVE_RANGES = "YES",
+                               GDAL_NUM_THREADS = "ALL_CPUS"
+                             )) {
+  args <- mget(names(formals()))
+  args$`...` <- NULL
+  args <- c(args, rlang::list2(...))
+  do.call(get_stac_data, args)
+}
+
 #' @rdname get_stac_data
 #' @export
 get_dem <- function(aoi,
diff --git a/man/get_stac_data.Rd b/man/get_stac_data.Rd
index 672b894..7a10fc2 100644
--- a/man/get_stac_data.Rd
+++ b/man/get_stac_data.Rd
@@ -5,6 +5,7 @@
 \alias{get_sentinel1_imagery}
 \alias{get_sentinel2_imagery}
 \alias{get_landsat_imagery}
+\alias{get_naip_imagery}
 \alias{get_dem}
 \title{Retrieve raster data from STAC endpoints}
 \usage{
@@ -117,6 +118,30 @@ get_landsat_imagery(
     GDAL_HTTP_MERGE_CONSECUTIVE_RANGES = "YES", GDAL_NUM_THREADS = "ALL_CPUS")
 )
 
+get_naip_imagery(
+  aoi,
+  start_date,
+  end_date,
+  ...,
+  pixel_x_size = 1,
+  pixel_y_size = 1,
+  asset_names = "image",
+  stac_source = "https://planetarycomputer.microsoft.com/api/stac/v1",
+  collection = "naip",
+  query_function = default_query_function,
+  sign_function = sign_planetary_computer,
+  rescale_bands = FALSE,
+  output_filename = paste0(proceduralnames::make_english_names(1), ".tif"),
+  composite_function = "merge",
+  limit = 999,
+  gdalwarp_options = c("-r", "bilinear", "-multi", "-overwrite", "-co",
+    "COMPRESS=DEFLATE", "-co", "PREDICTOR=2", "-co", "NUM_THREADS=ALL_CPUS"),
+  gdal_config_options = c(VSI_CACHE = "TRUE", GDAL_CACHEMAX = "30\%", VSI_CACHE_SIZE =
+    "10000000", GDAL_HTTP_MULTIPLEX = "YES", GDAL_INGESTED_BYTES_AT_OPEN = "32000",
+    GDAL_DISABLE_READDIR_ON_OPEN = "EMPTY_DIR", GDAL_HTTP_VERSION = "2",
+    GDAL_HTTP_MERGE_CONSECUTIVE_RANGES = "YES", GDAL_NUM_THREADS = "ALL_CPUS")
+)
+
 get_dem(
   aoi,
   ...,
diff --git a/tests/testthat/test-get_stac_data.R b/tests/testthat/test-get_stac_data.R
index 46a9f9d..887bdca 100644
--- a/tests/testthat/test-get_stac_data.R
+++ b/tests/testthat/test-get_stac_data.R
@@ -302,3 +302,21 @@ test_that("no-composites return the same data", {
     setdiff(rsi::landsat_band_mapping$planetary_computer_v1, "T")
   )
 })
+
+test_that("get_naip_imagery() is stable", {
+  skip_on_cran()
+  skip_if_offline()
+  aoi <- sf::st_point(c(-74.912131, 44.080410))
+  aoi <- sf::st_set_crs(sf::st_sfc(aoi), 4326)
+  aoi <- sf::st_buffer(sf::st_transform(aoi, 3857), 100)
+
+  expect_no_error(
+    out <- get_naip_imagery(
+      aoi,
+      "2018-01-01",
+      "2020-01-31",
+      output_filename = tempfile(fileext = ".tif")
+    )
+  )
+  expect_no_error(terra::rast(out))
+})