How to remove lakes/water from maps in R

Question

(I'll apologize in advance for not having a reprex. Can't get to work-going to be a separate question later today).

I have pulled ZCTA level data from the US Census, rolled the ZCTAs into groups, and created a choropleth map. I would like to remove the various lake boundaries. In places where the lake features are a larger portion of the total area (or are near the boundaries of my regions) their boundaries are a bit of an eyesore to the viewer.

library(tigris)
library(sf)
library(dplyr)
library(tidycensus)
library(stringr)
library(ggplot2)




var <- c(EduTotal = "B16010_001")
zip_sf <- get_acs(geography = "zcta",
                  variables = var,
                  year = 2017, survey = "acs5",
                  output = "wide", geometry = TRUE,
                  keep_geo_vars=TRUE
)

zip_sf %>%
  filter(str_detect(ZCTA5CE10,'^1')) %>%
    mutate(zip2=str_sub(ZCTA5CE10,1,2)) %>%
       group_by(zip2)  %>%
         summarize(meanEd=mean(EduTotalE))  %>%
  ggplot(aes(fill = meanEd)) + 
  geom_sf(col='red')

Zipcodes starting with '1'

One way to make this a reprex is to recognize that we don't need the actual census data, which requires having an API key, in order to work on the map, since that's what your problem actually is. But how do you plan on removing lakes? They're basically holes in the shapes—fill them in? — camille, Mar 13 '19 at 18:26
(Not sure I follow the comment on reprex. reprex is not working on any code I try it with-I just get errors) How to remove the lakes is my question. Not sure there is a way to do it, but there are smarter more experienced folks here that may have some ideas. — rclivornese, Mar 14 '19 at 16:10
Are your "lake boundaries" any rings internal to the external ring or rings defining your features? — Spacedman, Mar 17 '19 at 23:31

score 2 · Answer 1 · answered Mar 18 '19 at 00:07

What I think you want to do is draw the feature but only outline the external ring, or in the case of a feature with islands, all the external rings.

In PostGIS Simple Features theres an ST_ExternalRing function, but this doesn't seem implemented in sf package yet. You could ask Edzer nicely...

Meanwhile, this seems to work. Convert MULTIPOLYGON geometries to LINESTRINGS, convert those LINESTRINGS to POLYGONS and then UNION the polygons. In the process the holes (lakes) lose their identity as holes and the UNION process will drop them.

Example:

Run example(st_multipolygon) to create mp1 object. This is a MULTIPOLYGON object of three squares, two of which have holes in them:

> plot(mp1)

To drop the holes, do:

> mp1ext = 
    st_union(
      st_cast(
         st_cast(
            st_boundary(st_sfc(mp1)),
         "LINESTRING"),
      "POLYGON")
     )

Then to do your map, plot the area with a colour and a missing outline colour using the original data with the holes:

> plot(mp1,col="green",border=NA)

then add the borders without the holes using the new object:

> plot(mp1ext, col=NA, lwd=3,add=TRUE)

Note how the square holes (lakes) aren't outlined.

This is obviously an outline solution for a single object and uses base graphics instead of ggplot to do the plotting but the principles are probably here enough for you to adapt to your data. Other complications may arise. Write your own st_external_ring ring function based on my code and tweak it until it works better!

I think you are changing ring direction with your approach. If you use `st_sfc(st_multipolygon(lapply(mp1, "[", 1)))` to get the exterior rings, direction is preserved AFAIKT. — TimSalabim, Mar 18 '19 at 07:59

How to remove lakes/water from maps in R

1 Answers1