How to remove areas of water from the US census tract map using R/RStudio

Question

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_tracts_sans_water <- get_tracts("MI",year=2018)
mi_tracts_sans_water %>%
  ggplot(aes(fill = STATEFP)) +
  geom_sf()

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?

Answer 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")

---

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)

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)

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.