Merge pull request #96 from DOI-USGS/normalize_weights

for #94 -- normalized weights
DOI-USGS · Feb 3, 2024 · 56fda75 · 56fda75
2 parents 6a75800 + bf3da0a
commit 56fda75
Show file tree

Hide file tree

Showing 27 changed files with 633 additions and 67 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,7 +1,7 @@
 Package: ncdfgeom
 Type: Package
 Title: 'NetCDF' Geometry and Time Series
-Version: 1.1.6
+Version: 1.2.0
 Authors@R: c(person("David", "Blodgett", role = c("aut", "cre"), 
                     email = "dblodgett@usgs.gov"),
              person("Luke", "Winslow", role = "ctb"))

diff --git a/NAMESPACE b/NAMESPACE
@@ -21,6 +21,7 @@ importFrom(RNetCDF,var.inq.nc)
 importFrom(RNetCDF,var.put.nc)
 importFrom(dplyr,filter)
 importFrom(dplyr,group_by)
+importFrom(dplyr,left_join)
 importFrom(dplyr,mutate)
 importFrom(dplyr,right_join)
 importFrom(dplyr,select)
@@ -39,6 +40,7 @@ importFrom(sf,st_contains)
 importFrom(sf,st_convex_hull)
 importFrom(sf,st_coordinates)
 importFrom(sf,st_crs)
+importFrom(sf,st_drop_geometry)
 importFrom(sf,st_geometry)
 importFrom(sf,st_geometry_type)
 importFrom(sf,st_intersection)

diff --git a/NEWS.md b/NEWS.md
@@ -1,3 +1,7 @@
+1.2.0
+==========
+* calculate_area_intersection_weights() now supports normalized weights
+
 1.1.6
 ==========
 * remove geoknife as suggested package

diff --git a/R/calculate_area_intersection_weights.R b/R/calculate_area_intersection_weights.R
@@ -11,12 +11,22 @@
 #'
 #' @param x sf data.frame source features including one geometry column and one identifier column
 #' @param y sf data.frame target features including one geometry column and one identifier column
+#' @param normalize logical return normalized weights or not. See details and examples.
 #' @param allow_lonlat boolean If FALSE (the default) lon/lat target features are not allowed.
 #' Intersections in lon/lat are generally not valid and problematic at the international date line.
 #' @return data.frame containing fraction of each feature in x that is
-#' covered by each feature in y. e.g. If a feature from x is entirely within a feature in y,
-#' w will be 1. If a feature from x is 50% in one feature for y and 50% in another, there
-#' will be two rows, one for each x/y pair of features with w = 0.5 in each.
+#' covered by each feature in y. 
+#' 
+#' @details
+#' 
+#' Two versions of weights are available: 
+#' 
+#' `normalize = FALSE`, if a polygon from x is entirely within a polygon in y,
+#' w will be 1. If a polygon from x is 50% in one polygon from y and 50% in another, there
+#' will be two rows, one for each x/y pair of features with w = 0.5 in each. Weights
+#' will sum to 1 per SOURCE polygon if the target polygons fully cover that feature.
+#' `normalize = TRUE`, weights are divided by the target polygon area such that weights
+#' sum to 1 per TARGET polygon if the target polygon is fully covered by source polygons.
 #'
 #' @examples
 #' b1 = sf::st_polygon(list(rbind(c(-1,-1), c(1,-1),
@@ -40,14 +50,95 @@
 #'
 #' sf::st_agr(a) <- sf::st_agr(b) <- "constant"
 #'
-#' a_b <- calculate_area_intersection_weights(a, b)
-#' b_a <- calculate_area_intersection_weights(b, a)
+#' calculate_area_intersection_weights(a, b, normalize = FALSE)
+#' calculate_area_intersection_weights(a, b, normalize = TRUE)
+#' calculate_area_intersection_weights(b, a, normalize = FALSE)
+#' calculate_area_intersection_weights(b, a, normalize = TRUE)
+#' 
+#' #a more typical arrangement of polygons
+#' 
+#' b1 = sf::st_polygon(list(rbind(c(-1,-1), c(1,-1),
+#'                            c(1,1), c(-1,1),
+#'                            c(-1,-1))))
+#' b2 = b1 + 2
+#' b3 = b1 + c(0, 2)
+#' b4 = b1 + c(2, 0)
+#' b = sf::st_sfc(b1, b2, b3, b4)
+#' a1 = b1 * 0.75 + c(-.25, -.25)
+#' a2 = a1 + 1.5
+#' a3 = a1 + c(0, 1.5)
+#' a4 = a1 + c(1.5, 0)
+#' a = sf::st_sfc(a1,a2,a3, a4)
+#' plot(b, border = 'red', lwd = 3)
+#' plot(a, border = 'green', add = TRUE)
+#' 
+#' sf::st_crs(b) <- sf::st_crs(a) <- sf::st_crs(5070)
+#'
+#' b <- sf::st_sf(b, data.frame(idb = c(1, 2, 3, 4)))
+#' a <- sf::st_sf(a, data.frame(ida = c("a", "b", "c", "d")))
+#'
+#' sf::st_agr(a) <- sf::st_agr(b) <- "constant"
+#'
+#' # say we have data from `a` that we want sampled to `b`.
+#' # this gives the percent of each `a` that intersects each `b`
+#' 
+#' (a_b <- calculate_area_intersection_weights(a, b, normalize = FALSE))
+#' 
+#' # note that `w` sums to 1 where `b` completely covers `a`.
+#' 
+#' dplyr::summarize(dplyr::group_by(a_b, ida), w = sum(w))
+#' 
+#' # We can apply these weights like...
+#' dplyr::tibble(ida = unique(a_b$ida), 
+#'                    val = c(1, 2, 3, 4)) |>
+#'   dplyr::left_join(a_b, by = "ida") |>
+#'   dplyr::mutate(val = ifelse(is.na(w), NA, val),
+#'                 areasqkm = 1.5 ^ 2) |> # area of each polygon in `a`
+#'   dplyr::group_by(idb) |> # group so we get one row per `b`
+#'   # now we weight by the percent of the area from each polygon in `b` per polygon in `a`
+#'   dplyr::summarize(new_val = sum( (val * w * areasqkm), na.rm = TRUE ) / sum(w * areasqkm))
+#'   
+#' # we can go in reverse if we had data from b that we want sampled to a
+#' 
+#' (b_a <- calculate_area_intersection_weights(b, a, normalize = FALSE))
+#' 
+#' # note that `w` sums to 1 only where `a` complete covers `b`
+#' 
+#' dplyr::summarize(dplyr::group_by(b_a, idb), w = sum(w))
+#' 
+#' # with `normalize = TRUE`, `w` will sum to 1 when the target 
+#' # completely covers the source rather than when the source completely
+#' # covers the target. 
+#' 
+#' (a_b <- calculate_area_intersection_weights(a, b, normalize = TRUE))
+#' 
+#' dplyr::summarize(dplyr::group_by(a_b, idb), w = sum(w))
+#' 
+#' (b_a <- calculate_area_intersection_weights(b, a, normalize = TRUE))
+#' 
+#' dplyr::summarize(dplyr::group_by(b_a, ida), w = sum(w))
+#' 
+#' # We can apply these weights like...
+#' # Note that we don't need area in the normalized case
+#' dplyr::tibble(ida = unique(a_b$ida), 
+#'                    val = c(1, 2, 3, 4)) |>
+#'   dplyr::left_join(a_b, by = "ida") |>
+#'   dplyr::mutate(val = ifelse(is.na(w), NA, val)) |>
+#'   dplyr::group_by(idb) |> # group so we get one row per `b`
+#'   # now we weight by the percent of the area from each polygon in `b` per polygon in `a`
+#'   dplyr::summarize(new_val = sum( (val * w), na.rm = TRUE ))
+
 #'
 #' @export
-#' @importFrom sf st_intersection st_set_geometry st_area st_crs
-#' @importFrom dplyr mutate group_by right_join select ungroup
+#' @importFrom sf st_intersection st_set_geometry st_area st_crs st_drop_geometry
+#' @importFrom dplyr mutate group_by right_join select ungroup left_join mutate
 
-calculate_area_intersection_weights <- function(x, y, allow_lonlat = FALSE) {
+calculate_area_intersection_weights <- function(x, y, normalize, allow_lonlat = FALSE) {
+
+  if(missing(normalize)) {
+    warning("Required input normalize is missing, defaulting to FALSE.")
+    normalize <- FALSE
+  }
 
   if(!requireNamespace("areal")) stop("areal package required for intersection weights")
 
@@ -80,20 +171,39 @@ calculate_area_intersection_weights <- function(x, y, allow_lonlat = FALSE) {
 
   int <- areal::aw_intersect(y,
                              source = x,
-                             areaVar = "area")
+                             areaVar = "area_intersection")
   int <- areal::aw_total(int, source = x,
-                         id = "varx",
-                         areaVar = "area",
-                         totalVar = "totalArea",
+                         id = "varx", # the unique id in the "source" x
+                         areaVar = "area_intersection",
+                         totalVar = "totalArea_x",
                          type = "extensive",
                          weight = "total")
-  int <- areal::aw_weight(int, areaVar = "area",
-                          totalVar = "totalArea",
-                          areaWeight = "areaWeight")
+  int <- areal::aw_weight(int, areaVar = "area_intersection",
+                          totalVar = "totalArea_x",
+                          areaWeight = "areaWeight_x_y")
 
-  int <- right_join(st_set_geometry(int, NULL), st_set_geometry(x, NULL), by = "varx")
+  int <- right_join(st_drop_geometry(int), st_drop_geometry(x), by = "varx")
 
-  int <- select(int, varx, vary, w = areaWeight)
+  if(normalize) {
+
+    # for normalized, we return the percent of each target covered by each source
+    int <- areal::aw_intersect(y,
+                               source = x,
+                               areaVar = "area_intersection")
+
+    # for normalized, we sum the intersection area by the total target intersection area
+    int <- ungroup(mutate(group_by(int, vary), totalArea_y = sum(area_intersection)))
+
+    int <- areal::aw_weight(int, 
+                            areaVar = "area_intersection",
+                            totalVar = "totalArea_y",
+                            areaWeight = "areaWeight_x_y")
+
+  }
+
+  int <- right_join(st_drop_geometry(int), st_drop_geometry(x), by = "varx")
+
+  int <- select(int, varx, vary, w = areaWeight_x_y)
 
   names(int) <- c(id_x, id_y, "w")
 

diff --git a/docs/404.html b/docs/404.html
diff --git a/docs/DISCLAIMER.html b/docs/DISCLAIMER.html
diff --git a/docs/articles/geometry.html b/docs/articles/geometry.html
diff --git a/docs/articles/index.html b/docs/articles/index.html
diff --git a/docs/articles/ncdfgeom.html b/docs/articles/ncdfgeom.html
diff --git a/docs/articles/timeseries.html b/docs/articles/timeseries.html
diff --git a/docs/authors.html b/docs/authors.html
diff --git a/docs/index.html b/docs/index.html
diff --git a/docs/news/index.html b/docs/news/index.html
diff --git a/docs/pkgdown.yml b/docs/pkgdown.yml
@@ -5,5 +5,5 @@ articles:
   geometry: geometry.html
   ncdfgeom: ncdfgeom.html
   timeseries: timeseries.html
-last_built: 2024-01-22T18:39Z
+last_built: 2024-02-03T13:28Z
 
diff --git a/docs/reference/Rplot001.png b/docs/reference/Rplot001.png
diff --git a/docs/reference/calculate_area_intersection_weights-2.png b/docs/reference/calculate_area_intersection_weights-2.png