readme

Author: ramarty

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: blackmarbler
 Title: Black Marble Data and Statistics
-Version: 0.1.2
+Version: 0.2.0
 Authors@R: 
     c(person("Robert", "Marty", , "[email protected]", role = c("aut", "cre"),
            comment = c(ORCID = "0000-0002-3164-3813")),
@@ -20,7 +20,7 @@ Imports:
   purrr,
   lubridate,
   tidyr,
-  raster,
+  terra,
   sf,
   exactextractr,
   stringr,

R/blackmarbler.R

@@ -43,13 +43,6 @@ month_start_day_to_month <- function(x){
 
 month_start_day_to_month <- Vectorize(month_start_day_to_month)
 
-pad2 <- function(x){
-  if(nchar(x) == 1) out <- paste0("0", x)
-  if(nchar(x) == 2) out <- paste0(x)
-  return(out)
-}
-pad2 <- Vectorize(pad2)
-
 pad3 <- function(x){
   if(nchar(x) == 1) out <- paste0("00", x)
   if(nchar(x) == 2) out <- paste0("0", x)
@@ -313,12 +306,6 @@ file_to_raster <- function(f,
                       crs = myCrs,
                       extent = c(xMin,xMax,yMin,yMax))
 
-  #create extents class
-  #rasExt <- raster::extent(c(xMin,xMax,yMin,yMax))
-  
-  #assign the extents to the raster
-  #extent(outr) <- rasExt
-  
   #set fill values to NA
   outr <- remove_fill_value(outr, variable)
 
@@ -346,7 +333,7 @@ read_bm_csv <- function(year,
       df
     },
     error = function(e){
-      warning(paste0("Error with year: ", year, "; day: ", day))
+      #warning(paste0("Error with year: ", year, "; day: ", day))
       data.frame(NULL)
     }
   )
@@ -382,7 +369,7 @@ create_dataset_name_df <- function(product_id,
   }
 
   if(product_id == "VNP46A3"){
-    param_df <- cross_df(list(year            = 2012:year_end,
+    param_df <- cross_df(list(year = 2012:year_end,
                               day = c("001", "032", "061", "092", "122", "153", "183", "214", "245", "275", "306", "336",
                                       "060", "091", "121", "152", "182", "213", "244", "274", "305", "335")))
   }
@@ -422,10 +409,16 @@ create_dataset_name_df <- function(product_id,
   }
 
   #### Create data
-  files_df <- purrr::map2_dfr(param_df$year,
-                              param_df$day,
-                              read_bm_csv,
-                              product_id)
+  # files_df <- purrr::map2_dfr(param_df$year,
+  #                             param_df$day,
+  #                             read_bm_csv,
+  #                             product_id)
+  
+  files_df <- purrr::map2(param_df$year,
+                          param_df$day,
+                          read_bm_csv,
+                          product_id) %>%
+    bind_rows()
 
   return(files_df)
 }
@@ -533,7 +526,7 @@ count_n_obs <- function(values, coverage_fraction) {
 #' * `"VNP46A2"`: Daily (corrected)
 #' * `"VNP46A3"`: Monthly
 #' * `"VNP46A4"`: Annual
-#' @param date Date of raster data. Entering one date will produce a raster. Entering multiple dates will produce a raster stack.
+#' @param date Date of raster data. Entering one date will produce a `SpatRaster` object. Entering multiple dates will produce a `SpatRaster` object with multiple bands; one band per date.
 #' * For `product_id`s `"VNP46A1"` and `"VNP46A2"`, a date (eg, `"2021-10-03"`).
 #' * For `product_id` `"VNP46A3"`, a date or year-month (e.g., `"2021-10-01"`, where the day will be ignored, or `"2021-10"`).
 #' * For `product_id` `"VNP46A4"`, year or date  (e.g., `"2021-10-01"`, where the month and day will be ignored, or `2021`).
@@ -559,15 +552,15 @@ count_n_obs <- function(values, coverage_fraction) {
 #' - `1`: Poor-quality, The number of observations used for the composite is less than or equal to 3
 #' - `2`: Gap filled NTL based on historical data
 #' @param check_all_tiles_exist Check whether all Black Marble nighttime light tiles exist for the region of interest. Sometimes not all tiles are available, so the full region of interest may not be covered. If `TRUE`, skips cases where not all tiles are available. (Default: `TRUE`).
-#' @param interpol_na When data for more than one date is downloaded, whether to interpolate `NA` values in rasters using the `raster::approxNA` function. Additional arguments for the `raster::approxNA` function can also be passed into `bm_extract` (eg, `method`, `rule`, `f`, `ties`, `z`, `NA_rule`). (Default: `FALSE`).
+#' @param interpol_na When data for more than one date is downloaded, whether to interpolate `NA` values in rasters using the `terra::approximate` function. Additional arguments for the `terra::approximate` function can also be passed into `bm_extract` (eg, `method`, `rule`, `f`, `ties`, `z`, `NA_rule`). (Default: `FALSE`).
 #' @param output_location_type Where to produce output; either `memory` or `file`. If `memory`, functions returns a dataframe in R. If `file`, function exports a `.csv` file and returns `NULL`.
 #' @param file_dir (If `output_location_type = file`). The directory where data should be exported (default: `NULL`, so the working directory will be used)
 #' @param file_prefix (If `output_location_type = file`). Prefix to add to the file to be saved. The file will be saved as the following: `[file_prefix][product_id]_t[date].csv`
 #' @param file_skip_if_exists (If `output_location_type = file`). Whether the function should first check wither the file already exists, and to skip downloading or extracting data if the data for that date if the file already exists (default: `TRUE`).
 #' @param h5_dir Black Marble data are originally downloaded as `h5` files. If `h5_dir = NULL`, the function downloads to a temporary directory then deletes the directory. If `h5_dir` is set to a path, `h5` files are saved to that directory and not deleted. The function will then check if the needed `h5` file already exists in the directory; if it exists, the function will not re-download the `h5` file.
 #' @param quiet Suppress output that show downloading progress and other messages. (Default: `FALSE`).
 #'
-#' @param ... Additional arguments for `raster::approxNA`, if `interpol_na = TRUE`
+#' @param ... Additional arguments for `terra::approximate`, if `interpol_na = TRUE`
 #'
 #' @return Raster
 #'
@@ -628,6 +621,11 @@ bm_extract <- function(roi_sf,
     warning("interpol_na ignored. Interpolation only occurs when output_location_type = 'memory'")
   }
 
+  if(class(roi_sf)[1] == "SpatVector") roi_sf <- roi_sf %>% st_as_sf()
+  if(!("sf" %in% class(roi_sf))){
+    stop("roi must be an sf object")
+  }
+  
   # Assign interpolation variables ---------------------------------------------
   if(interpol_na == T){
     if(!exists("method")) method <- "linear"
@@ -676,13 +674,13 @@ bm_extract <- function(roi_sf,
                       quiet = quiet,
                       temp_dir = temp_dir)
 
-    bm_r <- raster::approxNA(bm_r,
-                             method = method,
-                             rule   = rule,
-                             f      = f,
-                             ties   = ties,
-                             z      = z,
-                             NArule = NArule)
+    bm_r <- terra::approximate(bm_r,
+                               method = method,
+                               rule   = rule,
+                               f      = f,
+                               ties   = ties,
+                               z      = z,
+                               NArule = NArule)
 
     #### Extract
     roi_df <- roi_sf %>% st_drop_geometry()
@@ -871,7 +869,7 @@ bm_extract <- function(roi_sf,
 #' * `"VNP46A2"`: Daily (corrected)
 #' * `"VNP46A3"`: Monthly
 #' * `"VNP46A4"`: Annual
-#' @param date Date of raster data. Entering one date will produce a raster. Entering multiple dates will produce a raster stack.
+#' @param date Date of raster data. Entering one date will produce a `SpatRaster` object. Entering multiple dates will produce a `SpatRaster` object with multiple bands; one band per date.
 #' * For `product_id`s `"VNP46A1"` and `"VNP46A2"`, a date (eg, `"2021-10-03"`).
 #' * For `product_id` `"VNP46A3"`, a date or year-month (e.g., `"2021-10-01"`, where the day will be ignored, or `"2021-10"`).
 #' * For `product_id` `"VNP46A4"`, year or date  (e.g., `"2021-10-01"`, where the month and day will be ignored, or `2021`).
@@ -895,15 +893,15 @@ bm_extract <- function(roi_sf,
 #' - `1`: Poor-quality, The number of observations used for the composite is less than or equal to 3
 #' - `2`: Gap filled NTL based on historical data
 #' @param check_all_tiles_exist Check whether all Black Marble nighttime light tiles exist for the region of interest. Sometimes not all tiles are available, so the full region of interest may not be covered. If `TRUE`, skips cases where not all tiles are available. (Default: `TRUE`).
-#' @param interpol_na When data for more than one date is downloaded, whether to interpolate `NA` values using the `raster::approxNA` function. Additional arguments for the `raster::approxNA` function can also be passed into `bm_raster` (eg, `method`, `rule`, `f`, `ties`, `z`, `NA_rule`). (Default: `FALSE`).
+#' @param interpol_na When data for more than one date is downloaded, whether to interpolate `NA` values using the `terra::approximate` function. Additional arguments for the `terra::approximate` function can also be passed into `bm_raster` (eg, `method`, `rule`, `f`, `ties`, `z`, `NA_rule`). (Default: `FALSE`).
 #' @param output_location_type Where to produce output; either `memory` or `file`. If `memory`, functions returns a raster in R. If `file`, function exports a `.tif` file and returns `NULL`.
 #' For `output_location_type = file`:
 #' @param file_dir The directory where data should be exported (default: `NULL`, so the working directory will be used)
 #' @param file_prefix Prefix to add to the file to be saved. The file will be saved as the following: `[file_prefix][product_id]_t[date].tif`
 #' @param file_skip_if_exists Whether the function should first check wither the file already exists, and to skip downloading or extracting data if the data for that date if the file already exists (default: `TRUE`).
 #' @param h5_dir Black Marble data are originally downloaded as `h5` files. If `h5_dir = NULL`, the function downloads to a temporary directory then deletes the directory. If `h5_dir` is set to a path, `h5` files are saved to that directory and not deleted. The function will then check if the needed `h5` file already exists in the directory; if it exists, the function will not re-download the `h5` file.
 #' @param quiet Suppress output that show downloading progress and other messages. (Default: `FALSE`).
-#' @param ... Additional arguments for `raster::approxNA`, if `interpol_na = TRUE`
+#' @param ... Additional arguments for `terra::approximate`, if `interpol_na = TRUE`
 #'
 #' @return Raster
 #'
@@ -1053,7 +1051,7 @@ bm_raster <- function(roi_sf,
             writeRaster(r, out_path)
 
           } else{
-            warning(paste0('"', out_path, '" already exists; skipping.\n'))
+            message(paste0('"', out_path, '" already exists; skipping.\n'))
           }
 
           r_out <- NULL # Saving as tif file, so output from function should be NULL
@@ -1090,20 +1088,20 @@ bm_raster <- function(roi_sf,
   if(length(r_list) == 1){
     r <- r_list[[1]]
   } else if (length(r_list) > 1){
-    r <- raster::stack(r_list)
+    r <- terra::rast(r_list)
   } else{
     r <- NULL
   }
 
   # Interpolate ----------------------------------------------------------------
   if(interpol_na %in% T){
-    r <- raster::approxNA(r,
-                          method = method,
-                          rule   = rule,
-                          f      = f,
-                          ties   = ties,
-                          z      = z,
-                          NArule = NArule)
+    r <- terra::approximate(r,
+                            method = method,
+                            rule   = rule,
+                            f      = f,
+                            ties   = ties,
+                            z      = z,
+                            NArule = NArule)
   }
 
   unlink(temp_dir, recursive = T)

README.md

@@ -18,7 +18,7 @@
 * [Usage](#usage)
   * [Setup](#setup)
   * [Make raster](#raster)
-  * [Make raster stack across multiple time periods](#stack)
+  * [Make raster across multiple time periods](#stack)
   * [Make map](#map)
   * [Make figure of trends in nighttime lights](#trends)
   * [Workflow to update data](#update-data)
@@ -71,8 +71,9 @@ Before downloading and extracting Black Marble data, we first load packages, def
 library(blackmarbler)
 library(geodata)
 library(sf)
-library(raster)
+library(terra)
 library(ggplot2)
+library(tidyterra)
 
 #### Define NASA bearer token
 bearer <- "BEARER-TOKEN-HERE"
@@ -81,7 +82,7 @@ bearer <- "BEARER-TOKEN-HERE"
 # Define region of interest (roi). The roi must be (1) an sf polygon and (2)
 # in the WGS84 (epsg:4326) coordinate reference system. Here, we use the
 # getData function to load a polygon of Ghana
-roi_sf <- gadm(country = "GHA", level=1, path = tempdir()) |> st_as_sf()
+roi_sf <- gadm(country = "GHA", level=1, path = tempdir()) 
 ```
 
 ### Make raster of nighttime lights <a name="raster">
@@ -109,9 +110,9 @@ r_2021 <- bm_raster(roi_sf = roi_sf,
                     bearer = bearer)
 ```
 
-### Make raster stack of nighttime lights across multiple time periods <a name="stack">
+### Make raster of nighttime lights across multiple time periods <a name="stack">
 
-To extract data for multiple time periods, add multiple time periods to `date`. The function will return a raster stack, where each raster band corresponds to a different date. The below code provides examples getting data across multiple days, months, and years.
+To extract data for multiple time periods, add multiple time periods to `date`. The function will return a `SpatRaster` object with multiple bands, where each band corresponds to a different date. The below code provides examples getting data across multiple days, months, and years.
 
 ```r
 #### Daily data in March 2021
@@ -145,28 +146,20 @@ r <- bm_raster(roi_sf = roi_sf,
                bearer = bearer)
 
 #### Prep data
-r <- r |> mask(roi_sf)
-
-r_df <- rasterToPoints(r, spatial = TRUE) |> as.data.frame()
-names(r_df) <- c("value", "x", "y")
-
-## Remove very low values of NTL; can be considered noise
-r_df$value[r_df$value <= 2] <- 0
+r <- r |> terraLLmask(roi_sf)
 
 ## Distribution is skewed, so log
-r_df$value_adj <- log(r_df$value+1)
+r[] <- log(r[] + 1)
 
 ##### Map
 p <- ggplot() +
-  geom_raster(data = r_df,
-  aes(x = x, y = y,
-  fill = value_adj)) +
+  geom_spatraster(data = r) +
   scale_fill_gradient2(low = "black",
                        mid = "yellow",
                        high = "red",
                        midpoint = 4.5) +
   labs(title = "Nighttime Lights: October 2021") +
-  coord_quickmap() +
+  coord_sf() +
   theme_void() +
   theme(plot.title = element_text(face = "bold", hjust = 0.5),
   legend.position = "none")
@@ -260,7 +253,7 @@ Both functions take the following arguments:
   * `"VNP46A3"`: Monthly
   * `"VNP46A4"`: Annual
 
-* **date:**  Date of raster data. Entering one date will produce a raster. Entering multiple dates will produce a raster stack.
+* **date:**  Date of raster data. Entering one date will produce a `SpatRaster` object. Entering multiple dates will produce a `SpatRaster` object with multiple bands; one band per date.
 
   * For `product_id`s `"VNP46A1"` and `"VNP46A2"`, a date (eg, `"2021-10-03"`).
   * For `product_id` `"VNP46A3"`, a date or year-month (e.g., `"2021-10-01"`, where the day will be ignored, or `"2021-10"`).
@@ -289,7 +282,7 @@ Both functions take the following arguments:
     * `2`: Gap filled NTL based on historical data
 
 * **check_all_tiles_exist:** Check whether all Black Marble nighttime light tiles exist for the region of interest. Sometimes not all tiles are available, so the full region of interest may not be covered. If `TRUE`, skips cases where not all tiles are available. (Default: `TRUE`).
-* **interpol_na:** When data for more than one date is downloaded, whether to interpolate `NA` values in rasters using the [`raster::approxNA`](https://www.rdocumentation.org/packages/raster/versions/3.6-26/topics/approxNA) function. Additional arguments for the [`raster::approxNA`](https://www.rdocumentation.org/packages/raster/versions/3.6-26/topics/approxNA) function can also be passed into `bm_raster`/`bm_extract` (eg, `method`, `rule`, `f`, `ties`, `z`, `NA_rule`). (Default: `FALSE`).
+* **interpol_na:** When data for more than one date is downloaded, whether to interpolate `NA` values in rasters using the [`terra::approximate`](https://www.rdocumentation.org/packages/raster/versions/3.6-26/topics/approxNA) function. Additional arguments for the [`terra::approximate`](https://www.rdocumentation.org/packages/raster/versions/3.6-26/topics/approxNA) function can also be passed into `bm_raster`/`bm_extract` (eg, `method`, `rule`, `f`, `ties`, `z`, `NA_rule`). (Default: `FALSE`).
 * **h5_dir:** Black Marble data are originally downloaded as `h5` files. If `h5_dir = NULL`, the function downloads to a temporary directory then deletes the directory. If `h5_dir` is set to a path, `h5` files are saved to that directory and not deleted. The function will then check if the needed `h5` file already exists in the directory; if it exists, the function will not re-download the `h5` file.
 
 
@@ -304,7 +297,7 @@ If `output_location_type = "file"`, the following arguments can be used:
 * **file_prefix:** Prefix to add to the file to be saved. The file will be saved as the following: `[file_prefix][product_id]_t[date].[tif/Rds]`
 * **file_skip_if_exists:** Whether the function should first check wither the file already exists, and to skip downloading or extracting data if the data for that date if the file already exists (default: `TRUE`). If the function is first run with `date = c(2018, 2019, 2020)`, then is later run with `date = c(2018, 2019, 2020, 2021)`, the function will only download/extract data for 2021. Skipping existing files can facilitate re-running the function at a later date to download only more recent data.
 
-* **...:** Additional arguments for [`raster::approxNA`](https://www.rdocumentation.org/packages/raster/versions/3.6-26/topics/approxNA), if `interpol_na = TRUE`
+* **...:** Additional arguments for [`terra::approximate`](https://rspatial.github.io/terra/reference/approximate.html), if `interpol_na = TRUE`
 
 ### Argument for `bm_extract` only <a name="args-extract">