chris-prener/areal

R package for areal interpolation

I am data scientist passionate about building real world evidence data and software products that accelerate the production of insights.

Currently, I work as a Director in Pfizer’s Evidence Generation Platform, where I am a RWE Digital Platforms specialist working across all of our theraputic areas.

I specialize in health disparities research, geospatial mapping and analyses, the development of open source research software, and data acquisition.

Follow My Work

• You can follow my work on my website

areal's Issues

Grid Interpolation

I think interpolating to grids is a frequent enough use case to warrant its own function. Below is a reprex of what I'm imagining.

# Helper Function to make grid from shape
# Requires that data is projected and that units are known

vars <- c("B02001_001E","B02001_003E")

test <- tidycensus::get_acs("tract", vars, output = "wide", state = "MO", county = 510, geometry = TRUE) %>%
  dplyr::transmute(black_pop = B02001_003E/B02001_001E,
            tid = row_number()) %>%
  sf::st_transform(102696)

make_grid <- function(.data){

  # find units of object
  units <- st_crs(.data, parameters = TRUE)$units_gdal

  # error if not projected
  if(units == "degree"){
    stop("Data must be projected to interpolate to grid squares")
  }

  cellsize <- switch (units,
    "Foot_US" = c(3280.84,3280.84), # Feet to Sq KM
    "Meter" = c(1000, 1000)
  ) # Will return NULL if no matches, need to implement more possible GDAL Units

  # make grid and clip to full boundary of data
  grid <- sf::st_make_grid(.data, cellsize, square = TRUE) %>%
    sf::st_intersection(sf::st_union(.data))

  # coerce to sf/df and add unique ID
  grid %<>% sf::st_sf() %>% dplyr::mutate(gid = dplyr::row_number())

  return(grid)
}

grid <- make_grid(test)
out <- areal::aw_interpolate(grid, "gid", test, "tid", weight = "sum", output = "sf", intensive = "black_pop")

Definitely needs refinement to implement, and the reason aw_interpolate is external to the function is because I could not get NSE to play nicely from within the function. Very neat too is that you can change from square to hexagons by setting square = FALSE

Pycnophylactic interpolation method

Adding functionality for this approach is a long term goal, and this issue will serve as a venue for discussion / documenting progress on it.

CRAN Release Checklist

Before we release to CRAN:

• Add examples to all exported functions
• Last round of spell checking
• Run package through devtools::release_checks() - passed
• Run package through r-hub - r-hub not used because it lacks dependencies needed to build sf and lwgeom on Debian
• Run package through Winbuilder
• Update NEWS.md
• Release what we've got a v0.1.2

Final steps:

• Create branch for CRAN release
• Update install sections on README.Rmd, index.Rmd, and vignettes/areal.Rmd
• Final R CMD check
• submit to CRAN!

Suggested approach for estimating interpolated margins of error with census data?

This may be outside the scope of the package, but I'm using aw_interpolate() to create demographic estimates for geographies that are not published by the Census/American Community Survey (e.g. community districts in New York City). Is there a suggested approach for aggregating the margins of error for the interpolated estimates?

When using entire census tracts or other geographies to derive new estimates, the Census Bureau provides recommended formulas for calculating the derived margins of error (implemented nicely in tidycensus). Is there an equivalent that could be used when interpolating?

weight and output arguments of aw_interpolate not picking up defaults?

Describe the bug
When weight and output arguments to aw_interpolate() are omitted, an error is returned.

Expected behavior
I assumed I could leave these blank and it would pick up the default settings in the function definition, which appear to be weight = "sum" and output = "sf". Not sure if this is the intended behavoir!

To Reproduce

areal::aw_interpolate(
  ar_stl_wards,
  tid = WARD,
  source = ar_stl_race,
  sid = GEOID,
  extensive = "TOTAL_E"
  )
#> Error in areal::aw_interpolate(ar_stl_wards, tid = WARD, source = ar_stl_race, : A weight type (either 'sum' or 'total') must be specified for the 'weight' argument.

^{Created on 2020-12-08 by the reprex package (v0.3.0)}

aw_interpolate fails when `tid` and `sid` are the same

See Issue #6 for a reprex - since @mfherman is using two census geometries, the ID is GEOID in both cases and we get this for output:

> aw_interpolate(
+   si_zcta,
+   tid = GEOID,
+   source = si_tract,
+   sid = GEOID,
+   weight = "sum",
+   output = "tibble",
+   extensive = "popE"
+ )
# A tibble: 13 x 2
   GEOID  popE
 * <chr> <dbl>
 1 10301    NA
 2 10302    NA
 3 10303    NA
 4 10304    NA
 5 10305    NA
 6 10306    NA
 7 10307    NA
 8 10308    NA
 9 10309    NA
10 10310    NA
11 10311    NA
12 10312    NA
13 10314    NA

Renaming one of the ID variables from GEOID to anything else fixes it. Not sure how to handle this yet. For now, @mfherman, rename your target ID variable until we address this.

Change to new cran checks badge URL

👋🏽 I maintain the cran checks badges. Please change to the new cran checks badge URL (e.g., https://badges.cranchecks.info/worst/dplyr.svg). Old badges at (e.g. https://cranchecks.info/badges/worst/dplyr) will be unavailable as of Jan 1st 2023.

• Blog post from June 2022 announcing the change: https://recology.info/2022/06/cran-checks-badges/
• New badges repo: https://github.com/r-hub/cchecksbadges
• Old badges repo: https://github.com/sckott/cchecksapi

Test Data

Need to add smaller test data for unit testing. And move full st louis data to a separate repository [To be accessed from gateway::gw_get_data()]

Test data should include:

Dasymetric

• A subset of Tracts containing a normal distribution of building density
• A corresponding subset of building footprints
• A corresponding target geometry

3-Class

• A subset of Tracts containing a normal distribution of LD and HD zoning
• A corresponding subset of zoning data
• A corresponding target geometry

Grammatical suggestions for the paper

Just a couple grammatical suggestions for the paper in
openjournals/joss-reviews#1221. Please adopt what you think fits and leave those that don't.

Introduction

Despite this, researchers often use them as proxies for neighborhoods since they come with readily available data and the techniques for estimating population values for overlapping but incongruent areas, such as a shapefile of neighborhood boundaries for a city, are less accessible.

This statement is fairly long, and it isn't entirely clear to me what the "are less accessible" is referring to. Would you mind trying to clarify this?

Unlike interpolation techniques, like inverse distance weighted interpolation, that are appropriate for point data, areal weighted interpolation is designed specifically for working with already aggregated data to some set of polygon spatial features

"like" is fairly informal for academic writing. Does "such as" or similar work?

The estimation process vary based on whether data are spatially extensive

Should be "varies" (or "estimation processes vary")

a census tract with a large park in it, for example, or a neighborhood that has a significant commercial area alongside residential housing

The grammar is a bit off here. How about: "for example, a census tract with a large park or a neighborhood that has a significant commercial area alongside residential housing would violate this assumption."

The areal R package

We also aim to provide additional functionality in comparison to the one existing approach found in the sf package [@pebesma2018simple], the st_interpolate_aw() function, while providing support for modern data management (e.g. tidyverse) and spatial data (e.g. sf) frameworks.

I think this can be lightly re-worded to be more clear: e.g.
"We also aim to provide additional functionality that is not currently available in existing approaches found in the sf package [@pebesma2018simple], for example the st_interpolate_aw() function, while providing support for modern data management (e.g. tidyverse) and spatial data (e.g. sf) frameworks."

allow for friction-less, iterative interpolation of both spatially extensive (e.g., count) and intensive (e.g., ratio) data, and

"frictionless" (no dash is needed)

`_pkgdown.yml` issue

The first line in the _pkgdown.yml file needs to have the url updated

3-class regression dasymetric method

Adding functionality for this approach is a long term goal, and this issue will serve as a venue for discussion / documenting progress on it.

aw_area assumes geometry column is named "geometry"

See here: https://github.com/slu-openGIS/areal/blob/master/R/aw_intersect.R#L103

The use of quosures is beyond my tidyverse skill or I'd PR a fix. You need to use sf::st_geometry() to operate on the geometry list column directly or lookup the name with attr(x, "sf_column").

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(geom = b, data.frame(idb = c(1, 2, 3, 4)))
a <- sf::st_sf(geom = a, data.frame(ida = c(1, 2, 3)))

sf::st_agr(a) <- sf::st_agr(b) <- "constant"

aw_intersect(a, b, "area")

attr(a, "sf_column") <- "geometry"
attr(b, "sf_column") <- "geometry"

names(a) <- c("ida", "geometry")
names(b) <- c("idb", "geometry")

aw_intersect(a, b, "area")

Upcoming sf release breaks areal

We see the following on R CMD check areal with upcoming sf 0.9-5:

  Running ‘testthat.R’
 ERROR
Running the tests in ‘tests/testthat.R’ failed.
Last 13 lines of output:
  Backtrace:
    1. testthat::expect_error(...)
    6. areal::aw_interpolate(...)
    7. areal:::aw_interpolate_single(...)
    8. areal:::aw_interpolater(...)
   13. areal::aw_intersect(., source = source, areaVar = "...area")
   18. areal:::aw_area(intersection, areaVar = !!areaVarQ)
   20. sf:::rename.sf(.data, geometry = colName)
  
  ══ testthat results  ═══════════════════════════════════════════════════════════
  [ OK: 145 | SKIPPED: 0 | WARNINGS: 0 | FAILED: 1 ]
  1. Failure: geom column is correctly renamed during interpolation (@test_aw_interpolate.R#151) 
  
  Error: testthat unit tests failed
  Execution halted
* checking for unstated dependencies in vignettes ... OK
* checking package vignettes in ‘inst/doc’ ... OK
* checking running R code from vignettes ...
  ‘areal-weighted-interpolation.Rmd’ using ‘UTF-8’... OK
  ‘areal.Rmd’ using ‘UTF-8’... OK
  ‘data-preparation.Rmd’ using ‘UTF-8’... OK
 NONE
* checking re-building of vignette outputs ... WARNING
Error(s) in re-building vignettes:
  ...
--- re-building ‘areal-weighted-interpolation.Rmd’ using rmarkdown
--- finished re-building ‘areal-weighted-interpolation.Rmd’

--- re-building ‘areal.Rmd’ using rmarkdown
--- finished re-building ‘areal.Rmd’

--- re-building ‘data-preparation.Rmd’ using rmarkdown
Quitting from lines 177-180 (data-preparation.Rmd) 
Error: processing vignette 'data-preparation.Rmd' failed with diagnostics:
'names' attribute [5] must be the same length as the vector [4]
--- failed re-building ‘data-preparation.Rmd’

SUMMARY: processing the following file failed:
  ‘data-preparation.Rmd’

Error: Vignette re-building failed.
Execution halted

* checking PDF version of manual ... OK
* DONE

Status: 1 ERROR, 1 WARNING
See
  ‘/tmp/areal.Rcheck/00check.log’
for details.

which is probably due to r-spatial/sf#1431 - please anticipate fixing & updating (or considering not checking at this granularity, as these may be issues areal should not care about)

Transition to Github Actions

Our CI Pipeline needs to get upgraded as we look to integrate all of these new functions. Can be nearly carbon-copied from censusxy

I'll do this whenever we discuss.

aw_interpolate() fails when attempting to cast to MULTIPOLYGON

Describe the bug
aw_interpolate() fails when when attempting to cast geometry to MULTIPOLYGON. Possibly due to source and target sf objects being MULTIPOLYGON and not POLYGON?

To Reproduce
I'm not sure if this is the most minimal example as it involves grabbing data and geometry using tidycensus, but I believe it reproduced the error I was having working with the other shapefiles I was working with. I first try aw_interpolate() and then break down that function to pinpoint where it seems like the error is occurring for me.

library(tidyverse)
library(tidycensus)
library(areal)
library(sf)
#> Linking to GEOS 3.7.0, GDAL 2.3.2, PROJ 5.2.0

si_tract <- get_acs(
  geography = "tract",
  variables = c("pop" = "B01001_001"),
  state = "NY",
  county = "Richmond",
  geometry = TRUE,
  output = "wide"
  ) %>% 
  st_transform(2263) %>% 
  select(GEOID, popE)
#> Getting data from the 2012-2016 5-year ACS

si_zcta <- get_acs(
  geography = "zip code tabulation area",
  variables = c("pop" = "B01001_001"),
  geometry = TRUE,
  output = "wide"
  ) %>% 
  filter(GEOID %in% 10301:10314) %>% 
  st_transform(2263) %>% 
  select(GEOID)
#> Getting data from the 2012-2016 5-year ACS

ar_validate(si_tract, si_zcta, varList = "popE")
#> [1] TRUE

aw_interpolate(
  si_zcta,
  tid = GEOID,
  source = si_tract,
  sid = GEOID,
  weight = "sum",
  output = "tibble",
  extensive = "popE"
)
#> Warning in st_cast.GEOMETRYCOLLECTION(X[[i]], ...): only first part of
#> geometrycollection is retained
#> Error in st_cast.POINT(x[[1]], to, ...): cannot create MULTIPOLYGON from POINT

aw_intersect(si_zcta, si_tract, areaVar = "area")
#> Warning in st_cast.GEOMETRYCOLLECTION(X[[i]], ...): only first part of
#> geometrycollection is retained
#> Error in st_cast.POINT(x[[1]], to, ...): cannot create MULTIPOLYGON from POINT

intersection <- st_intersection(si_tract, si_zcta)
#> Warning: attribute variables are assumed to be spatially constant
#> throughout all geometries
out <- areal:::aw_area(intersection, "area")
st_cast(out, "MULTIPOLYGON")
#> Warning in st_cast.GEOMETRYCOLLECTION(X[[i]], ...): only first part of
#> geometrycollection is retained
#> Error in st_cast.POINT(x[[1]], to, ...): cannot create MULTIPOLYGON from POINT

^{Created on 2018-12-20 by the reprex package (v0.2.1)}

Desktop (please complete the following information):

• OS: macOS 10.14.1
• Version of R 3.5.1
• Version of RStudio 1.2.1194

Additional context
I'm able to successful run the examples in the vignettes on this machine.

Implement Inverse Distance Weighting Algorithm

Leaving this as a note to self when I have more time to get back to areal.

Inverse distance weighting is an invaluable tool for creating raster surfaces from point samples. However, the implementations of this technique in R are poorly documented and rely on (highly problematic) Spatial* class objects. I've written, but not yet perfected, an sf implementation after spending hours trying to replicate sp methods.

It seems appropriate that an sf-based abstraction of IDW gets added as a function to areal

Tibble Output Missing Other Variables

When output = "tibble" is selected, areal fails to return other variables in the original data.

aw_interpolate/ar_validate fail despite planar CRS

I'm trying to do a simple areal interpolation of data coming straight from tidycensus and tigris, both transformed to crs = 26915, but despite these transformations the interpolation fails, and using ar_validate show "CRS is Planar" as "FALSE". Code shown below, any idea why this is failing?

blocks.shp <- transform(tigris::blocks(state = state, county = counties,year = 2021,class = "sf"),crs = 26915)
bgs.shp <- transform(tigris::block_groups(state = state, county = counties,year = 2021,class = "sf"),crs = 26915) %>% left_join(y = se,by = "GEOID")
se_shp <- aw_interpolate(blocks.shp,tid = "GEOID20",source = bgs.shp,sid = "GEOID", extensive = c("TotalPop","MinPop","NMinPop","LIPop","NLIPop","LI185"), output = "sf")

Missing packages in vignette

Need to add a couple of packages to the vignette to reduce new user confusion - see this RStudio Community post

Interpolating Average (Gaussian) Data

Hey @chris-prener! This package is awesome, and has been incredibly helpful for working with census data at the neighborhood level. Currently I've been using it and exploring some analyses with count data - but many variables of census data that I'd like to explore are Gaussian data (median income, average year built).

Do you have any plans to include a way to interpolate these types of data in the areal package?
Currently, as also described in your package documentation, the workflow for conducting an interpolation for data types like this is extensive and not easily reproducible in ArcGIS.

I also admittedly know that I may be asking for a feature that may not be possible because Gaussian data makes a certain set of assumptions that may not be as easily put into this package - I also could be making this more complex and this type of data would be explored using the intensive option.

Oddity with aw_intersects and adding up areas?

Running a few tests with aw_intersect, aw_total, and aw_weight, I noticed that if I run the three in a pipeline for a single target feature and a large feature coverage of source data, I get back weights that are all 1. The reason for this appears to be that the intersect step is cutting off all the source features outside of the target polygon as shown in this screenshot.

As you can see here: https://github.com/dblodgett-usgs/intersectr/blob/master/R/calculate_area_intersection_weights.R#L70 I use the area of the source polygons themselves to get that area...

Not sure which is right here? But worth noting anyways. Or am I missing something that needs to be documented or I missed in the documentation?

GitHub Actions

Tasks:

1. Set-up CI build on GitHub actions
2. Address build errors
3. Add pkgdown GitHub action for site build
4. README.md and repo housekeeping (remove old CI files, update badges)
5. Ensure code coverage is working correctly

Add polygon checks to ar_validate()

While using this dev branch I wasted an inordinate amount of time having not checked my intermediate shapefiles geometries that weren't POLYGON or MULTIPOLYGON.

Would adding the following check add too much overhead to aw_dasymetric?

if(any(!st_geometry_type(nyc_buildings_tidy) %in% c("POLYGON", "MULTIPOLYGON"))){
  stop('The `intermediate` shapefiles must contain only `POLYGON` and `MULTIPOLYGON` geometries. Remove other geometries before passing to aw_dasymetric')
}

Originally posted by @charliejhadley in #27 (comment)

Namespaces in Imports field not imported from

When I run R CMD check on the package, I see

* checking dependencies in R code ... NOTE
Namespaces in Imports field not imported from:
  ‘lwgeom’ ‘tibble’

suggesting that these two packages can be removed from Imports:
openjournals/joss-reviews#1221

install error

@chris-prener when installing via

> install.packages("areal")

I'm getting:

ERROR: dependencies ‘lwgeom’, ‘sf’ are not available for package ‘areal’
  * removing ‘/home/serge/anaconda3/lib/R/library/areal’

The downloaded source packages are in
        ‘/tmp/RtmpBEpBpJ/downloaded_packages’
Updating HTML index of packages in '.Library'
Making 'packages.html' ... done
Warning messages:
1: In install.packages("areal") :
  installation of package ‘units’ had non-zero exit status
2: In install.packages("areal") :
  installation of package ‘sf’ had non-zero exit status
3: In install.packages("areal") :
  installation of package ‘lwgeom’ had non-zero exit status
4: In install.packages("areal") :
  installation of package ‘areal’ had non-zero exit status
> library(areal)
Error in library(areal) : there is no package called ‘areal’
>

I realize this might be due to using R installed via conda. But on Ubuntu 18.04 I'm getting the same error(s) on a clean install of R from the official repositories:

* DONE (dplyr)
ERROR: dependencies ‘lwgeom’, ‘sf’ are not available for package ‘areal’
* removing ‘/home/serge/R/x86_64-pc-linux-gnu-library/3.4/areal’

The downloaded source packages are in
        ‘/tmp/Rtmpb2kPO5/downloaded_packages’
Warning messages:
1: In install.packages("areal") :
  installation of package ‘units’ had non-zero exit status
2: In install.packages("areal") :
  installation of package ‘sf’ had non-zero exit status
3: In install.packages("areal") :
  installation of package ‘lwgeom’ had non-zero exit status
4: In install.packages("areal") :
  installation of package ‘areal’ had non-zero exit status
> library(areal)
Error in library(areal) : there is no package called ‘areal’

Intelligent Dasymetric Mapping

I would like to be able to reproduce the results of this toolbox in R:

https://github.com/USEPA/Dasymetric-Toolbox-ArcGISPro,

Binary dasymetric method

Adding functionality for this approach is a long term goal, and this issue will serve as a venue for discussion / documenting progress on it.

dplyr 1.0.0

We're seeing these errors when testing areal against the soon released dplyr 1.0.0:

[master*] 126 MiB ❯ revdepcheck::revdep_details(, "areal")
══ Reverse dependency check ═════════════════════════════════════ areal 0.1.5 ══

Status: BROKEN

── Newly failing

✖ checking tests ...

── Before ──────────────────────────────────────────────────────────────────────
0 errors ✔ | 0 warnings ✔ | 0 notes ✔

── After ───────────────────────────────────────────────────────────────────────
❯ checking tests ...
  See below...

── Test failures ───────────────────────────────────────────────── testthat ────

> library(testthat)
> library(areal)
>
> test_check("areal")
── 1. Failure: correctly specified functions execute without error (@test_aw_int
`aw_interpolate(...)` threw an error.
Message: Join columns must be present in data
✖ Problem with `WARD`
Class:   rlang_error/error/condition
Backtrace:
  1. testthat::expect_error(...)
  6. areal::aw_interpolate(...)
 12. dplyr:::left_join.data.frame(.data, data, by = tidQN)
 13. dplyr:::join_mutate(...)
 14. dplyr:::join_cols(...)
 15. dplyr:::standardise_join_by(by, x_names = x_names, y_names = y_names)
 16. dplyr:::check_join_vars(by$y, y_names)

── 2. Failure: correctly specified functions execute without error (@test_aw_int
`aw_interpolate(...)` threw an error.
Message: Join columns must be present in data
✖ Problem with `WARD`
Class:   rlang_error/error/condition
Backtrace:
  1. testthat::expect_error(...)
  6. areal::aw_interpolate(...)
 12. dplyr:::left_join.data.frame(.data, data, by = tidQN)
 13. dplyr:::join_mutate(...)
 14. dplyr:::join_cols(...)
 15. dplyr:::standardise_join_by(by, x_names = x_names, y_names = y_names)
 16. dplyr:::check_join_vars(by$y, y_names)

── 3. Failure: correctly specified functions execute without error (@test_aw_int
`aw_interpolate(...)` threw an error.
Message: Join columns must be present in data
✖ Problem with `WARD`
Class:   rlang_error/error/condition
Backtrace:
  1. testthat::expect_error(...)
  6. areal::aw_interpolate(...)
  8. dplyr:::left_join.data.frame(exresults, inresults, by = tidQN)
  9. dplyr:::join_mutate(...)
 10. dplyr:::join_cols(...)
 11. dplyr:::standardise_join_by(by, x_names = x_names, y_names = y_names)
 12. dplyr:::check_join_vars(by$x, x_names)

══ testthat results  ═══════════════════════════════════════════════════════════
[ OK: 138 | SKIPPED: 0 | WARNINGS: 6 | FAILED: 3 ]
1. Failure: correctly specified functions execute without error (@test_aw_interpolate.R#134)
2. Failure: correctly specified functions execute without error (@test_aw_interpolate.R#136)
3. Failure: correctly specified functions execute without error (@test_aw_interpolate.R#138)

Error: testthat unit tests failed
Execution halted

1 error ✖ | 0 warnings ✔ | 0 notes ✔

I have not fully investigated yet what's happening, but perhaps you'd know what is happening ? I can help if you would narrow it down a bit ...