Keep polygons that are fully or partially within a larger one without a negative buffer

user_12211397 0
Viewed 325

I'd like to keep the census tracts that comprise a larger census geography. The tracts should fit within the larger geography perfectly.

I can do this with a negative buffer, but is there a better way to do this?

library(tidyverse)
library(sf)
library(tidycensus)
library(tigris)
library(tmap)

balt <- places("Maryland",
               year = 2018,
               class = "sf") %>% 
  filter(NAME == "Baltimore") %>% 
  st_transform(4326)

balt_tracts <- get_acs(geography = "tract",
                       variables = "B01003_001E", # population
                       year = 2018,
                       state = "Maryland",
                       survey = "acs5",
                       geometry = TRUE) %>% 
  st_transform(4326)

balt_tracts %>% 
  .[st_within(st_buffer(., -0.001),  
              balt) %>% lengths > 0,] %>% 
  qtm()

enter image description here

An intersection returns tracts outside the larger geography that share an edge, which I don't want.

balt_tracts[balt,] %>% 
  qtm()

enter image description here

Edit

I'm looking for a way to take this result of st_overlaps() and keep only those polygons that are inside the boundary (balt).

balt_tracts[st_overlaps(balt_tracts, balt) %>% lengths > 0,] %>% 
  tm_shape() + tm_polygons(col = "GEOID")

enter image description here

One thing I did notice is that st_touches() returns only three polygons. Is this a precision issue, or does this mean the census tract borders are not actually perfectly aligned with the larger polygon's borders?

 balt_tracts[st_touches(balt_tracts, balt) %>% lengths > 0,] %>% 
   tm_shape() + tm_polygons(col = "GEOID")

enter image description here

sf r
2 Answers

In your specific example, you should just supply a county argument to get_acs(), as Baltimore city is an independent entity:

library(tidyverse)
library(sf)
library(tidycensus)
library(tigris)
library(tmap)

balt_tracts <- get_acs(geography = "tract",
                       variables = "B01003_001E", # population
                       year = 2018,
                       state = "Maryland",
                       county = "Baltimore city",
                       survey = "acs5",
                       geometry = TRUE) 

qtm(balt_tracts)

enter image description here

Your broader question is a worthwhile one, however. If you know that your geometries align perfectly (as they will when using tigris/tidycensus for the same year, provided that cb = TRUE or cb = FALSE is used consistently), you can use sf::st_filter() with the st_within spatial predicate to do spatial subsetting. For example, let's say we want to grab Census tracts in the Baltimore metro area (which is entirely in Maryland):

balt_metro <- core_based_statistical_areas(cb = TRUE, 
                                           year = 2018, 
                                           class = "sf") %>%
  filter(str_detect(NAME, "Baltimore"))

balt_metro_tracts <- tracts("MD", cb = TRUE, 
                            year = 2018, class = "sf") %>%
  st_filter(balt_metro, .predicate = st_within)

ggplot() + 
  geom_sf(data = balt_metro_tracts, fill = "white") + 
  geom_sf(data = balt_metro, fill = NA, color = "red")

enter image description here

user_56829620

One alternative is to turn each tract into a point and intersect on those, which removes the coincident borders that gave me trouble. This approach does not require choosing an arbitrary buffer distance like the one in the initial question.

st_point_on_surface() is preferable to st_centroid() since it guarantees the point will lie on the surface of the polygon, whereas some polygons will have centroids lying outside its surface.

balt_tracts[st_point_on_surface(balt_tracts)[balt,],] %>% 
  qtm()

enter image description here

user_122113970
Related