1

I am trying to aggregate all of the census tracts in the United States for analysis and display, including a realistic depiction of water features when zoomed in. I run into memory issues when doing this at scale, so I am looking for a more efficient method.

Here is my current function for doing this:

library(tidyverse)
library(tigris)
library(sf)
library(rmapshaper)

get_tracts <- function(state_abbr, 
                       year=NULL, 
                       refresh=FALSE, 
                       remove_water=TRUE){
  print(paste0("getting shapefile for tracts in ",state_abbr))
  state_tracts <- tracts(cb = TRUE, 
                         year = year, 
                         class="sf", 
                         state=state_abbr, 
                         refresh=refresh) %>% 
    st_transform(., crs=4326) %>% st_make_valid()
  Sys.sleep(1) #included for politeness to census servers
  if(!remove_water){
    return(state_tracts)
  }
  print(paste0("fetching water boundaries for ",state_abbr))

  county_codes <- unique(state_tracts$COUNTYFP)
  
  print(paste0("joining water areas into one shapefile for ",state_abbr))
  
  # https://community.rstudio.com/t/better-ways-to-remove-areas-of-water-from-us-map/60771
  # https://gis.stackexchange.com/questions/308119/function-in-r-to-get-difference-between-two-sets-of-polygons-comparable-in-speed
  water_area <- 
    map_dfr(
      county_codes,
      ~ area_water(state_abbr, 
                   county = ., 
                   year = year,
                   refresh = refresh,
                   class = "sf")  %>% 
        select("geometry") %>% 
        st_transform(., crs=4326) %>% 
        st_make_valid()
    )
  
  print("subtracting water areas from census tract areas")
  
  state_tracts_sans_water <- ms_erase(ms_simplify(state_tracts),
                                      ms_simplify(water_area), 
                                      remove_slivers=TRUE,
                                      sys=FALSE)
  Sys.sleep(1)
  return(state_tracts_sans_water)
}

This works fine for a state such as Michigan:

mi_tracts <- get_tracts("MI",year=2018,remove_water = FALSE)
mi_tracts %>%
  ggplot(aes(fill = STATEFP)) +
  geom_sf()

MI unfiltered census tracts

mi_tracts_sans_water <- get_tracts("MI",year=2018)
mi_tracts_sans_water %>%
  ggplot(aes(fill = STATEFP)) +
  geom_sf()

MI filtered to exclude bodies of water

However, for a large state like Texas (the state with the most counties), it fails by crashing my R session:

tx_tracts_sans_water <- get_tracts("TX",year=2018)

Sometimes it simply crashes the R session in RStudio, and when run as part of a loop of many states I may receive the following error:

V8 FATAL ERROR in Ineffective mark-compacts near heap limit: Allocation failed - JavaScript heap out of memory # # Fatal error in , line 0 # API fatal error handler returned after process out of memory # # # #FailureMessage Object: 0x7ffe249dd920 ==== C stack trace ===============================

Things I have thought about or tried:

  • Upgraded from standard sf functions to rmapshaper per this question/answer, which improved performance to the extent that I could even have the aforementioned issue.
  • Was unable to figure out how to use the system mapshaper using the sys=TRUE option of ms_erase on my RStudio Server setup, and question whether this would be sustainable given the need to define PATH specific to each user and version of node in use.
  • I am open to parallelization but due to the functions being performed I am not sure that this will be productive.
  • Batching the logic by county rather than by state. I am asking the question here before refactoring my code to do this, but if this is the only way or there are other pro's to this that make it preferable then I am open to it.
  • The answer to @mrblobby's RStudio Community post is to directly download the relevant hydrology shapefiles from the census website. However, this is not reproducible/automated. Automated ways to do this that can be parametrized by state are welcome.
  • To allow for zooming into specific communities at a high resolution, I would prefer not to filter out bodies of water based on their size (AWATER). However, if there is a standard that exists for such filtering then I am open to doing this.

So, with all of that being said, how should I go about removing areas of water from the US census tract map using R/RStudio?

1
  • 1
    FYI, not related to your problem but, in your function, correct R coding practice is to not use print but rather something like message("getting data for tracts in ", state_abbr). This also negates the use of paste. Also, when thinking about print consider cat as an alternative. This becomes quite relevant when submitting packages to CRAN. The package reviewer will catch inapproapriace function usage and request that you change it before the package is approved. ALso always invoke the NAMESPACE of functions used in your function (eg., tigris::tracts) Feb 22 at 15:32
1

The best solution is to include a basemap and increase the transparency of the data layer:

tx_tracts <- get_tracts("TX",year=2018,remove_water = FALSE)
library(ggmap)
# https://stackoverflow.com/questions/47749078/how-to-put-a-geom-sf-produced-map-on-top-of-a-ggmap-produced-raster/50844502#50844502
ggmap_bbox <- function(map) {
  if (!inherits(map, "ggmap")) stop("map must be a ggmap object")
  # Extract the bounding box (in lat/lon) from the ggmap to a numeric vector, 
  # and set the names to what sf::st_bbox expects:
  map_bbox <- setNames(unlist(attr(map, "bb")), 
                       c("ymin", "xmin", "ymax", "xmax"))

  # Coonvert the bbox to an sf polygon, transform it to 3857, 
  # and convert back to a bbox (convoluted, but it works)
  bbox_3857 <- st_bbox(st_transform(st_as_sfc(st_bbox(map_bbox, crs = 4326)), 3857))

  # Overwrite the bbox of the ggmap object with the transformed coordinates 
  attr(map, "bb")$ll.lat <- bbox_3857["ymin"]
  attr(map, "bb")$ll.lon <- bbox_3857["xmin"]
  attr(map, "bb")$ur.lat <- bbox_3857["ymax"]
  attr(map, "bb")$ur.lon <- bbox_3857["xmax"]
  map
}
b <- st_bbox(tx_tracts,crs=3857)
names(b) <- c("left","bottom","right","top")
# https://stackoverflow.com/questions/52704695/is-ggmap-broken-basic-qmap-produces-arguments-imply-differing-number-of-rows
tx_basemap <- get_stamenmap(bbox=b, 
                            zoom=calc_zoom(b, adjust=as.integer(0)),
                            maptype="toner-lite")#"terrain-background")

tx_basemap <- ggmap_bbox(tx_basemap)

ggmap(tx_basemap) + 
  geom_sf(data=st_transform(tx_tracts, 3857), 
          inherit.aes = FALSE, aes(fill=COUNTYFP), alpha=0.25, col=NA) +
  theme(legend.position = "none")

texas basemap example


Note that an overlay remains over the water, which might prompt a desire to still remove the water from the overlay. If you really want to remove the water then it is just a matter of doing so in small enough chunks as to not overburden the memory.

get_water <- function(state_abbr,
                  county,
                  year,
                  refresh,
                  class="sf"){
the_water <- area_water(state_abbr, 
             county = county, 
             year = year,
             refresh = refresh,
             class = "sf")  %>% 
    select("geometry") %>% 
    add_column(., COUNTYFP=county,.before=1) %>% 
    st_transform(., crs=4326)
  return(the_water)
}


get_tracts <- function(state_abbr, 
                       year=NULL, 
                       refresh=FALSE, 
                       remove_water=TRUE){
  # if we have saved them, then load them if refresh is FALSE
  print(paste0("getting shapefile for tracts in ",state_abbr))
  state_tracts <- tracts(cb = TRUE, 
                         year = year, 
                         class="sf", 
                         state=state_abbr, 
                         refresh=refresh) %>% 
    st_transform(., crs=4326) %>% st_make_valid()
  Sys.sleep(1)
  if(!remove_water){
    return(state_tracts)
  }
  print(paste0("fetching water boundaries for ",state_abbr))
  
  county_codes <- unique(state_tracts$COUNTYFP)
  
  print(paste0("joining water areas into one shapefile for ",state_abbr))
  
  # https://community.rstudio.com/t/better-ways-to-remove-areas-of-water-from-us-map/60771
  # https://gis.stackexchange.com/questions/308119/function-in-r-to-get-difference-between-two-sets-of-polygons-comparable-in-speed
  water_area <- 
    map_dfr(
      county_codes,
      ~get_water(
        state_abbr=state_abbr,
        county=.,
        year=year,
        refresh=refresh,
        class="sf"
      )
    )
  print(paste0("subtracting water from census tracts for ",state_abbr))
  pb <- txtProgressBar(min=0, 
                       max=nrow(state_tracts),
                       style=3)
  current_county_num <- -1
  for(i in 1:nrow(state_tracts)){
    next_county_num <- as.numeric(st_drop_geometry(state_tracts)[i,c("COUNTYFP")])
    if(next_county_num!=current_county_num){
      county_water_area <- water_area %>%
        filter(`COUNTYFP`==next_county_num)
      current_county_num <- next_county_num
    }
    if(nrow(county_water_area)>0){
      state_tracts[i,] <- ms_erase(state_tracts[i,], 
                                   county_water_area, 
                                   remove_slivers=TRUE)
    }
    setTxtProgressBar(pb, i)
  }
  return(state_tracts)
}

tx_tracts_sans_water <- get_tracts("TX",year=2018)
tx_tracts_sans_water %>%
  ggplot(aes(fill = STATEFP)) +
  geom_sf(col=NA)

texas cutout example

This takes a long time but it does "work." I did explore parallelizing this operation but did not have any luck with adapting the existing solutions (parallel, multidplyr) to spatial dataframes. With any luck, someone with expertise in one of those libraries will provide a better solution, or a data.table expert will show us how to do this in one line of code.

However, it's worth noting that this does not improve the aesthetic significantly because the basemap may not have the water visible at certain zoom levels. Or, the water being cutout may not quite match with the water as depicted on the basemap at certain zoom levels:

ggmap(tx_basemap) + 
  geom_sf(data=st_transform(tx_tracts_sans_water, 3857), 
          inherit.aes = FALSE, aes(fill=GEOID), alpha=0.5, col=NA)

texas cutout with basemap

Therefore, for my use case, I am simply using a basemap per the top of this answer, with the zoom level of the basemap defined by the extent of the geography into which I am zooming.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.