Merge pull request #48 from tcarleton/terra_conversion

Terra conversion

Author: cullen-molitor

.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/data_nj_counties.R

@@ -2,6 +2,31 @@
 # This is because we are referencing column names in dplyr / data.table objects which use
 # 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", "."))
+globalVariables(
+  c(
+    "cell_area_km2",
+    "coverage_fraction",
+    "day",
+    "era5_grid",
+    "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/global_variables.R

@@ -8,6 +8,7 @@
 #' @param polygon_id_col the name of a column in the `polygons` object with a
 #'   unique identifier for each polygon
 #' @param grid a raster layer with the same spatial resolution as the data
+#'   (must use geographic coordinates)
 #' @param secondary_weights an optional table of secondary weights, output from
 #'   the `secondary_weights()` function
 #'
@@ -16,7 +17,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 +33,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 +44,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
   ## -----------------------------------------------
@@ -62,6 +67,12 @@ overlay_weights <- function(polygons, polygon_id_col, grid = era5_grid, secondar
   rast_xmax <- terra::ext(clim_area_raster)$xmax
   rast_res <-  terra::xres(clim_area_raster)
 
+  ## check if SpatRaster is in geographic coodrinates
+  if(!terra::is.lonlat(clim_raster)) {
+    stop(crayon::red('Grid does not have geographic coordinates.'))
+
+  }
+
  ## stop if polygons are not in standard coordinate system
  if(poly_xmax > 180) {
 
@@ -70,7 +81,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(clim_raster) != 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 +108,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 +144,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)]
 
     }
 
@@ -165,25 +187,19 @@ overlay_weights <- function(polygons, polygon_id_col, grid = era5_grid, secondar
     w_merged[, weight := weight * w_area]
 
     # Create column that determines if entire polygon has a weight == 0
-    zero_polys <- data.frame(w_merged) |>
-      dplyr::group_by(poly_id) |>
-      dplyr::summarise(sum_weight = sum(weight)) |>
-      dplyr::ungroup() |>
-      dplyr::filter(sum_weight == 0) |>
-      dplyr::select(poly_id) |>
-      dplyr::distinct()
+    zero_polys <- w_merged[, .(sum_weight = sum(weight)),
+                                by = .(poly_id)]
+
+    zero_polys <- unique(zero_polys[sum_weight == 0, .(poly_id)])
 
-    if(nrow(zero_polys > 0)) {
+    if(nrow(zero_polys) > 0) {
 
       warning(crayon::red("Warning: weight = 0 for all pixels in some of your polygons; NAs will be returned for weights"))
 
     }
 
     # List any polygons with NA values in 1 or more grid cells
-    na_polys <- data.frame(w_merged) |>
-      dplyr::filter(is.na(weight)) |>
-      dplyr::select(poly_id) |>
-      dplyr::distinct()
+    na_polys <- unique(w_merged[is.na(weight), .(poly_id)])
 
     # # Warning if there are polygons with NA weight values
     # if(nrow(na_polys > 0)) {
@@ -193,12 +209,8 @@ overlay_weights <- function(polygons, polygon_id_col, grid = era5_grid, secondar
     # }
 
     # Update the weight to NA for all grid cells in na_polys
-    w_merged <- w_merged |>
-      dplyr::mutate(weight = ifelse(poly_id %in% c(na_polys$poly_id, zero_polys$poly_id), NA, weight)) |>
-      data.table::as.data.table()
-
-    ## this doesn't work with dt... figure out or delete and use above
-    # w_merged[, weight := data.table::fifelse(poly_id %in% c(na_polys$poly_id, zero_polys$poly_id), NA, weight)]
+    w_merged <- w_merged[, weight := ifelse(poly_id %in% c(na_polys$poly_id,
+                                                           zero_polys$poly_id), NA, weight)]
 
   }
 
@@ -231,7 +243,7 @@ overlay_weights <- function(polygons, polygon_id_col, grid = era5_grid, secondar
   # Check that polygon weights sum to 1 or 0 if all weights are NA
   if (!is.null(secondary_weights)){
 
-    for(i in nrow(check_weights)){
+    for(i in seq_len(nrow(check_weights))){
 
       if(!dplyr::near(check_weights$w_area[i], 1, tol=0.001)) {
 
@@ -246,7 +258,7 @@ overlay_weights <- function(polygons, polygon_id_col, grid = era5_grid, secondar
 
     } else {
 
-    for(i in nrow(check_weights)){
+    for(i in seq_len(nrow(check_weights))){
 
       if(!dplyr::near(check_weights$w_sum[i], 1, tol=0.001)){
 
@@ -262,10 +274,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/overlay_weights.R

@@ -7,7 +7,7 @@
 #' @param secondary_raster a raster of a secondary variable, for example
 #'   cropland coverage or population
 #' @param grid a raster layer with the same spatial resolution as the climate
-#'   data
+#'   data (must use geographic coordinates)
 #' @param extent an optional extent to crop the secondary_raster to for faster
 #'   processing. Format must be compatible with raster::crop(). The default is "full", which
 #'   resamples the whole secondary raster without cropping.
@@ -32,7 +32,7 @@ secondary_weights <- function(secondary_raster, grid = era5_grid, extent = "full
   ## Return error if terra::extent can't inherit from the value supplied
   ## won't be able to check if secondary raster fully overlaps if
   ## this input isn't compatible with terra::extent()
-  if(!is.character(extent) | length(extent) > 1){
+  if(!is.character(extent) || length(extent) > 1){
     tryCatch({
       terra::ext(extent)
       message(crayon::green("User-defined extent compatible with raster"))
@@ -72,12 +72,26 @@ 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
+
+  }
+
+  ## check if SpatRaster is in geographic coodrinates
+  if(!terra::is.lonlat(clim_raster)) {
+    stop(crayon::red('Grid does not have geographic coordinates.'))
+
+  }
+
   # 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 +151,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 %% 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 +180,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 %% 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.'))
 
@@ -186,14 +212,14 @@ secondary_weights <- function(secondary_raster, grid = era5_grid, extent = "full
   ## crop the ERA/climate raster to the appropriate extent
   ## use the extent of the previously user-cropped secondary raster
   ## -----------------------------------------------
+
   # Find the difference between the climate raster resolution and secondary raster resolution
   clim_raster <- terra::crop(clim_raster, terra::ext(secondary_raster), snap="out")
 
   ## set crs of secondary raster to match climate data
   ## -----------------------------------------------
   terra::crs(secondary_raster) <- terra::crs(clim_raster)
 
-
   ## check if the cropped secondary raster contains NA values
   if(isTRUE(any(is.na(terra::values(secondary_raster, na.rm=FALSE))))) {
 
@@ -206,7 +232,7 @@ secondary_weights <- function(secondary_raster, grid = era5_grid, extent = "full
   ## -----------------------------------------------
   message(crayon::green("Resampling secondary_raster"))
 
-  resampled_raster = terra::resample(secondary_raster, clim_raster, method="bilinear")
+  resampled_raster <- terra::resample(secondary_raster, clim_raster, method="bilinear")
 
   ## Make a data.table of the values of the resampled raster with lat/lon
   ## -----------------------------------------------

R/secondary_weights.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