Terra conversion

.Rbuildignore

@@ -5,3 +5,7 @@
 ^data-raw$
 ^doc$
 ^Meta$
+^_pkgdown\.yml$
+^docs$
+^pkgdown$
+^\.github$

.gitignore

@@ -6,3 +6,4 @@
 inst/doc
 /doc/
 /Meta/
+docs

DESCRIPTION

@@ -23,16 +23,19 @@ Imports:
     lubridate,
     magrittr,
     methods,
-    raster,
     sf,
     stringr,
     terra
 Suggests: 
+    climateR (>= 0.3.5),
     knitr,
+    raster,
     rmarkdown,
     testthat (>= 3.0.0),
     tidyr,
-    tigris,
-    usethis
+    usethis,
+    tigris
 Config/testthat/edition: 3
 VignetteBuilder: knitr
+Remotes:  
+    mikejohnson51/climateR

NAMESPACE

@@ -19,3 +19,4 @@ importFrom(magrittr,"%>%")
 importFrom(methods,as)
 importFrom(methods,setMethod)
 importFrom(methods,signature)
+importFrom(stats,na.omit)

R/data_cropland_world_2015_era5.R

@@ -1,6 +1,11 @@
 #' Global cropland weights
 #'
 #' A data.table with ERA5 resolution returned by running secondary_weights() with cropland data from
-#' 2015
+#' 2015 from https://glad.umd.edu/dataset/croplands
+#'
+#' Dataset Reference: P. Potapov, S. Turubanova, M.C. Hansen, A. Tyukavina, V.
+# Zalles, A. Khan, X.-P. Song, A. Pickens, Q. Shen, J. Cortez. (2021) Global
+# maps of cropland extent and change show accelerated cropland expansion in the
+# twenty-first century. Nature Food. https://doi.org/10.1038/s43016-021-00429-z
 #'
 "cropland_world_2015_era5"

R/global_variables.R

@@ -3,5 +3,7 @@
 # Non-Standard Evaluation.
 
 globalVariables(c("cell_area_km2", "coverage_fraction", "day", "era5_grid",
-                  "hour", "month", "poly_id", "sum_weight", "value", "w_area",
-                  "w_sum", "weight", "x", "y", "year", "."))
+                  "hour", "month", "minute", "poly_id", "sum_weight", "value",
+                  "w_area", "w_sum", "weight", "x", "y", "year", ".",
+                  "is_right_xmin", "is_left_xmax", "x_low", "x_high", "y_low",
+                  "y_high", "..cols_to_keep"))

R/overlay_weights.R

@@ -16,7 +16,7 @@
 #'
 #' @examples
 #' overlay_output_with_secondary_weights <- overlay_weights(
-#'   polygons = nj_counties, # Polygons outlining the 21 counties of New Jersey
+#'   polygons = tigris::counties("nj"), # Polygons outlining the 21 counties of New Jersey
 #'   polygon_id_col = "COUNTYFP", # The name of the column with the unique
 #'                                # county identifiers
 #'   grid = era5_grid, # The grid to use when extracting area weights (era5_grid is the
@@ -32,7 +32,7 @@
 #'
 #'
 #' overlay_output_without_secondary_weights <- overlay_weights(
-#'   polygons = nj_counties, # Polygons outlining the 21 counties of New Jersey
+#'   polygons = tigris::counties("nj"), # Polygons outlining the 21 counties of New Jersey
 #'   polygon_id_col = "COUNTYFP" # The name of the column with the unique county
 #'                               # identifiers
 #'   )
@@ -43,10 +43,14 @@
 #' @export
 overlay_weights <- function(polygons, polygon_id_col, grid = era5_grid, secondary_weights = NULL){
 
-  # Create raster (Putting terra::rast() here creates, for unknown reasons,
-  # issues with devtools::check())
-  clim_raster <- as(grid, "SpatRaster") # only reads the first band
+  ## check to make sure climate raster is a spatraster, change if not
+  if (!inherits(grid, "SpatRaster")) {
+    clim_raster <- terra::rast(grid)
+  } else {
+
+    clim_raster <- grid
 
+  }
 
   ## Raster cell area
   ## -----------------------------------------------
@@ -70,7 +74,13 @@ overlay_weights <- function(polygons, polygon_id_col, grid = era5_grid, secondar
  }
 
   ## check if coordinate systems match, if no shift raster to -180 to 180
-  if(rast_xmax > 180 + rast_res / 2) {
+  if(rast_xmax > 180 + rast_res) {
+
+    # Make sure the cell widths aren't peculiar otherwise the rotate function will
+    # mess things up
+    if(360 %% terra::xres(grid) != 0){
+      stop(crayon::red('Grid is in climate coordinate system (longitude 0 to 360) and grid cell width does not divide 360 evenly, making accurate alignment impossible.'))
+    }
 
     message(crayon::yellow('Aligning longitudes to standard coordinates.'))
 
@@ -91,8 +101,13 @@ overlay_weights <- function(polygons, polygon_id_col, grid = era5_grid, secondar
     ## rotate
     clim_area_raster <- terra::rotate(clim_area_raster)
 
-    }
+  }
 
+  # Extend the grid to cover all polygons consistent with the extended rotate
+  # above. exact_extract already does this in the background to a certain
+  # degree, so this just allows us to be explicit about how we handle NAs later
+  # on.
+  clim_area_raster <- terra::extend(clim_area_raster, terra::ext(polygons), snap = 'out')
 
   ## Match raster and polygon crs
   crs_raster <- terra::crs(clim_area_raster)
@@ -122,11 +137,11 @@ overlay_weights <- function(polygons, polygon_id_col, grid = era5_grid, secondar
     s_weight_max <- max(weights_dt$x)
 
     ## if secondary_weights is in 0-360, adjust x val
-    if(s_weight_max > 180 + rast_res) {
+    if(s_weight_max > 180 + rast_res / 2) {
 
       message(crayon::yellow('Adjusting secondary weights longitude to standard coordinates.'))
 
-      weights_dt[, x := data.table::fifelse(x > 180 + rast_res, x - 360, x)]
+      weights_dt[, x := data.table::fifelse(x > 180 + rast_res / 2, x - 360, x)]
 
     }
 
@@ -262,10 +277,15 @@ overlay_weights <- function(polygons, polygon_id_col, grid = era5_grid, secondar
   # If it doesn't error out then all weight sums = 1
   message(crayon::green('All weights sum to 1.'))
 
-  ## Convert back to 0-360
-  ## -----------------------------------------------
+  ## Return table in coordinate system that matches that of the climate data
+  ## ------------------------------------------------------------------------
+
+  if(rast_xmax > 180 + rast_res) {
+
+    w_norm[, x := data.table::fifelse(x < 0 + rast_res / 2, x + 360, x)]
+
+  }
 
-  w_norm[, x := data.table::fifelse(x < 0, x + 360, x)]
 
   return(w_norm)
 

R/secondary_weights.R

@@ -72,12 +72,20 @@ secondary_weights <- function(secondary_raster, grid = era5_grid, extent = "full
 
   }
 
+  ## check to make sure climate raster is a spatraster, change if not
+  if (!inherits(grid, "SpatRaster")) {
+    clim_raster <- terra::rast(grid)
+  } else {
+
+    clim_raster <- grid
+
+  }
+
   # Create climate raster from input raster
-  clim_raster <- terra::rast(grid) # only reads the first band
+  clim_raster <- clim_raster[[1]] # only reads the first band
 
   ## climate raster information for creating buffer and doing checks/rotations
   c_rast_xmax <- terra::ext(clim_raster)$xmax
-  # c_rast_xmin <- raster::extent(clim_raster)@xmin
 
   ## find xy resolution for rasters
   c_rast_xres <- terra::xres(clim_raster)
@@ -137,7 +145,13 @@ secondary_weights <- function(secondary_raster, grid = era5_grid, extent = "full
 
 
   ## if secondary raster in -180 to 180 and clim raster 0-360, rotate clim raster
-  if(s_rast_xmax <= (180 + s_rast_xres / 2) & c_rast_xmax >= (180 + c_rast_xres / 2)) {
+  if(s_rast_xmax < (180 + s_rast_xres) & c_rast_xmax > (180 + c_rast_xres)) {
+
+    # Make sure the cell widths aren't peculiar otherwise the rotate function will
+    # mess things up
+    if(360 %% terra::xres(c_rast_xres) != 0){
+      stop(crayon::red('Grid is in climate coordinate system (longitude 0 to 360) and grid cell width does not divide 360 evenly, making accurate alignment impossible.'))
+    }
 
     message(crayon::yellow('Longitude coordinates do not match. Aligning longitudes to standard coordinates.'))
 
@@ -160,7 +174,13 @@ secondary_weights <- function(secondary_raster, grid = era5_grid, extent = "full
   }
 
   ## if secondary raster in 0-360 and clim raster -180 to 180, rotate secondary raster
-  if(s_rast_xmax >= (180 + s_rast_xres / 2) & c_rast_xmax <= (180 + c_rast_xres / 2)) {
+  if(s_rast_xmax > (180 + s_rast_xres) & c_rast_xmax < (180 + c_rast_xres)) {
+
+    # Make sure the cell widths aren't peculiar otherwise the rotate function will
+    # mess things up
+    if(360 %% terra::xres(s_rast_xres) != 0){
+      stop(crayon::red('Grid is in climate coordinate system (longitude 0 to 360) and grid cell width does not divide 360 evenly, making accurate alignment impossible.'))
+    }
 
     message(crayon::yellow('Longitude coordinates do not match. Aligning longitudes to standard coordinates.'))
 

R/stagg-package.R

@@ -11,6 +11,7 @@
 #' @importFrom data.table .NGRP
 #' @importFrom data.table .SD
 #' @importFrom data.table data.table
+#' @importFrom stats na.omit
 #' @importFrom magrittr %>%
 #' @importFrom methods as
 #' @importFrom methods setMethod