diff --git a/DESCRIPTION b/DESCRIPTION
index 44e292c184d30e595190114c004182321f96bbf0..51ff5db59480603d0ce22c407baf585fe4952b8b 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -12,7 +12,7 @@ BugReports: https://github.com/USGS-R/ncdfgeom/issues
Imports: RNetCDF, ncmeta, sf, dplyr, methods, stars
Depends:
R (>= 3.5)
-Suggests: testthat, knitr, rmarkdown, pkgdown, geoknife, ncdf4, jsonlite
+Suggests: testthat, knitr, rmarkdown, pkgdown, geoknife, ncdf4, jsonlite, areal
License: CC0
Encoding: UTF-8
VignetteBuilder: knitr
diff --git a/NAMESPACE b/NAMESPACE
index 6ace2e3fbead8fb518b70a3df7c881a1ce650244..bf2aed6dd492cd8b8e0d984dbaf62f3e8b648581 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -1,5 +1,6 @@
# Generated by roxygen2: do not edit by hand
+export(calculate_area_intersection_weights)
export(create_cell_geometry)
export(read_attribute_data)
export(read_geometry)
@@ -20,13 +21,17 @@ importFrom(RNetCDF,var.inq.nc)
importFrom(RNetCDF,var.put.nc)
importFrom(dplyr,filter)
importFrom(dplyr,group_by)
+importFrom(dplyr,mutate)
+importFrom(dplyr,right_join)
importFrom(dplyr,select)
+importFrom(dplyr,ungroup)
importFrom(methods,is)
importFrom(ncmeta,nc_atts)
importFrom(ncmeta,nc_axes)
importFrom(ncmeta,nc_dim)
importFrom(ncmeta,nc_meta)
importFrom(sf,"st_crs<-")
+importFrom(sf,st_area)
importFrom(sf,st_as_sf)
importFrom(sf,st_bbox)
importFrom(sf,st_buffer)
@@ -36,6 +41,7 @@ importFrom(sf,st_coordinates)
importFrom(sf,st_crs)
importFrom(sf,st_geometry)
importFrom(sf,st_geometry_type)
+importFrom(sf,st_intersection)
importFrom(sf,st_join)
importFrom(sf,st_linestring)
importFrom(sf,st_multilinestring)
diff --git a/R/calculate_area_intersection_weights.R b/R/calculate_area_intersection_weights.R
new file mode 100644
index 0000000000000000000000000000000000000000..650d184c8a969500696a40f4c56479ef6f46036b
--- /dev/null
+++ b/R/calculate_area_intersection_weights.R
@@ -0,0 +1,103 @@
+#' Area Weighted Intersection (areal implementation)
+#' @description Returns the fractional percent of each
+#' feature in x that is covered by each intersecting feature
+#' in y. These can be used as the weights in an area-weighted
+#' mean overlay analysis where x is the data source and area-
+#' weighted means are being generated for the target, y.
+#'
+#' This function is a lightwieght wrapper around the functions
+#' \link[areal]{aw_intersect} \link[areal]{aw_total} and \link[areal]{aw_weight}
+#' from the \href{https://github.com/slu-openGIS/areal}{areal package}.
+#'
+#' @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 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.
+#'
+#' @examples
+#' 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, 2)
+#' b4 = b1 + c(2.2, 0)
+#' b = sf::st_sfc(b1, b2, b3, b4)
+#' a1 = b1 * 0.8
+#' a2 = a1 + c(1, 2)
+#' a3 = a1 + c(-1, 2)
+#' a = sf::st_sfc(a1,a2,a3)
+#' plot(b, border = 'red')
+#' 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(1, 2, 3)))
+#'
+#' 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)
+#'
+#' @export
+#' @importFrom sf st_intersection st_set_geometry st_area st_crs
+#' @importFrom dplyr mutate group_by right_join select ungroup
+
+calculate_area_intersection_weights <- function(x, y, allow_lonlat = FALSE) {
+
+ if(!requireNamespace("areal")) stop("areal package required for intersection weights")
+
+ if (st_crs(x) != st_crs(y)) {
+ x <- st_transform(x, st_crs(y))
+ }
+
+ if(st_crs(y)$proj == "longlat" & !allow_lonlat) {
+ stop("Found lon/lat coordinates and allow_lonlat is FALSE.")
+ }
+
+ # Standard evaluation is for chumps.
+ id_x <- names(x)[names(x) != attr(x, "sf_column")]
+ id_y <- names(y)[names(y) != attr(y, "sf_column")]
+
+ # There is a bug in areal and this works around it.
+ geom_name_x <- attr(x, "sf_column")
+ attr(x, "sf_column") <- "geometry"
+ names(x)[which(names(x) == geom_name_x)] <- "geometry"
+
+ geom_name_y <- attr(y, "sf_column")
+ attr(y, "sf_column") <- "geometry"
+ names(y)[which(names(y) == geom_name_y)] <- "geometry"
+
+ if (length(id_x) != 1 | length(id_y) != 1)
+ stop("x and y must have one and only one non-geometry column")
+
+ names(x)[names(x) == id_x] <- "varx"
+ names(y)[names(y) == id_y] <- "vary"
+
+ int <- areal::aw_intersect(y,
+ source = x,
+ areaVar = "area")
+ int <- areal::aw_total(int, source = x,
+ id = "varx",
+ areaVar = "area",
+ totalVar = "totalArea",
+ type = "extensive",
+ weight = "total")
+ int <- areal::aw_weight(int, areaVar = "area",
+ totalVar = "totalArea",
+ areaWeight = "areaWeight")
+
+ int <- right_join(st_set_geometry(int, NULL), st_set_geometry(x, NULL), by = "varx")
+
+ int <- select(int, varx, vary, w = areaWeight)
+
+ names(int) <- c(id_x, id_y, "w")
+
+ return(dplyr::as_tibble(int))
+}
+
+varx <- vary <- w <- d <- poly_id <- areaWeight <- NULL
diff --git a/man/calculate_area_intersection_weights.Rd b/man/calculate_area_intersection_weights.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..f93a1fa7049af8aba8915fc45f5b0254cfb32c28
--- /dev/null
+++ b/man/calculate_area_intersection_weights.Rd
@@ -0,0 +1,59 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/calculate_area_intersection_weights.R
+\name{calculate_area_intersection_weights}
+\alias{calculate_area_intersection_weights}
+\title{Area Weighted Intersection (areal implementation)}
+\usage{
+calculate_area_intersection_weights(x, y, allow_lonlat = FALSE)
+}
+\arguments{
+\item{x}{sf data.frame source features including one geometry column and one identifier column}
+
+\item{y}{sf data.frame target features including one geometry column and one identifier column}
+
+\item{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.}
+}
+\value{
+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.
+}
+\description{
+Returns the fractional percent of each
+feature in x that is covered by each intersecting feature
+in y. These can be used as the weights in an area-weighted
+mean overlay analysis where x is the data source and area-
+weighted means are being generated for the target, y.
+
+This function is a lightwieght wrapper around the functions
+\link[areal]{aw_intersect} \link[areal]{aw_total} and \link[areal]{aw_weight}
+from the \href{https://github.com/slu-openGIS/areal}{areal package}.
+}
+\examples{
+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, 2)
+b4 = b1 + c(2.2, 0)
+b = sf::st_sfc(b1, b2, b3, b4)
+a1 = b1 * 0.8
+a2 = a1 + c(1, 2)
+a3 = a1 + c(-1, 2)
+a = sf::st_sfc(a1,a2,a3)
+plot(b, border = 'red')
+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(1, 2, 3)))
+
+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)
+
+}
diff --git a/tests/testthat/test_area_weighted.R b/tests/testthat/test_area_weighted.R
new file mode 100644
index 0000000000000000000000000000000000000000..98c1a250ab6e79b84822218fa92cb238d8e6bb36
--- /dev/null
+++ b/tests/testthat/test_area_weighted.R
@@ -0,0 +1,49 @@
+context("calculate_area_intersection_weights")
+
+test_that("area_intersection", {
+
+ 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, 2)
+ b4 <- b1 + c(2.2, 0)
+ b <- sf::st_sfc(b1, b2, b3, b4)
+ a1 <- b1 * 0.8
+ a2 <- a1 + c(1, 2)
+ a3 <- a1 + c(-1, 2)
+ a <- sf::st_sfc(a1, a2, a3)
+
+ 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(1, 2, 3)))
+
+ 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)
+
+ expect_equal(as.numeric(a_b[1, ]),
+ c(1, 1, 1), info = "a1 is 100% covered by b1.")
+ expect_equal(as.numeric(a_b[2, ]),
+ c(2, 2, 0.5), info = "a2 is 50% covered by b2.")
+ expect_equal(as.numeric(a_b[3, ]),
+ c(2, 3, 0.375), info = "a2 is 37.5% covered by b3.")
+ expect_equal(as.numeric(a_b[4, ]),
+ c(3, 3, 0.625), info = "a3 is 62.5% covered by b3.")
+
+ expect_equal(as.numeric(b_a[1, ]),
+ c(1, 1, 0.64), info = "b1 is 64% covered by a1")
+ expect_equal(as.numeric(b_a[2, ]),
+ c(2, 2, 0.32), info = "b2 is 32% covered by a2")
+ expect_equal(as.numeric(b_a[3, ]),
+ c(3, 2, 0.24), info = "b3 is 24% covered by a2")
+ expect_equal(as.numeric(b_a[4, ]),
+ c(3, 3, 0.4), info = "b3 is 40% covered by a3")
+ expect_equal(data.frame(b_a[5, ]),
+ data.frame(tibble::tibble(idb = 4,
+ ida = as.numeric(NA),
+ w = as.numeric(NA))),
+ info = "b4 is not covered")
+})