aw_interpolate() fails when attempting to cast to MULTIPOLYGON

[mfherman]

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.

[mfherman]

Also, just wanted to give an early thanks for this package! I've done some areal interpolation in R before, but it makes a lot of sense to turn it into a package.

[chris-prener]

Thanks @mfherman!

This might be a bit of an edge case - the only reference I can find to the error from sf::st_cast(), which was the root case, is here.

Your intersection is creating geometry collections, which our test data has not done. Reading between the lines of the sf issue, it seems like you might be creating slivers with the intersection but I didn't verify that. I've added a line to aw_intersect() that gets a version of your code working correctly.

I'm going to leave this open until we merge the PR - we've got a backlog since the package is under CRAN review. What we'll do is wait until we hear back from the initial review. If we have to resubmit, will merge the multipolygon branch before we do.

Otherwise, will merge into master after CRAN acceptance and then submit a patch to CRAN in January.

[mfherman]

Thanks, @chris-prener! aw_intersect() on the multipolygon branch is working for my sf objects now.

[chris-prener]

Glad to hear it @mfherman!

[dblodgett-usgs]

@chris-prener -- I think the fix to this issue might be causing some performance issues when it isn't really needed.

A little back of the envelope profiling of my areal function vs a non-areal function gives me ~9 minutes for areal and ~30 seconds for non-areal. I dug into the function calls I'm making and found that this line is the killer. I also get the warning:

Warning message:
In st_collection_extract.sf(intersection) : x is already of type POLYGON.

Reproduce with: big_test.rda.zip

load("big_test.rda")
int <- areal::aw_intersect(target_polygons, data_source_cells, "area")

I'm not really sure what the consequences of not calling that are or what the best check would be to see if you don't need it... will do some research and get back to you.