0

I am using getData("GADM") to obtain map data for the US and Canada. Hudson's Bay and the Atlantic are transparent so I can add a blue background that shows through. However, all the Great Lakes have disappeared. Can I get them to be transparent as well to let the blue show through?

The data from GADM is detailed and the plotting of the maps is very slow in RStudio on my Mac. I have added code found here to speed up the plotting.

library(tidyverse)
library(raster)

# Convert the polygons into data frames so we can make lines
poly2df <- function(poly) {
  # Convert the polygons into data frames so we can make lines
  # Number of regions
  n_regions <- length(poly@polygons)

  # Get the coords into a data frame
  poly_df <- c()
  for(i in 1:n_regions) {
    # Number of polygons for first region
    n_poly <- length(poly@polygons[[i]]@Polygons)
    print(paste("There are",n_poly,"polygons"))
    # Create progress bar
    pb <- txtProgressBar(min = 0, max = n_poly, style = 3)
    for(j in 1:n_poly) {
      poly_df <- rbind(poly_df, NA, 
                       poly@polygons[[i]]@Polygons[[j]]@coords)
      # Update progress bar
      setTxtProgressBar(pb, j)
    }
    close(pb)
    print(paste("Finished region",i,"of",n_regions))
  }
  poly_df <- data.frame(poly_df)
  names(poly_df) <- c('lon','lat')
  return(poly_df)
}

# Get the main polygons, will determine by area.
getSmallPolys <- function(poly, minarea = 0.01) {
  # Get the areas
  areas <- lapply(poly@polygons, 
                  function(x) sapply(x@Polygons, function(y) y@area))

  # Quick summary of the areas
  print(quantile(unlist(areas)))

  # Which are the big polygons?
  bigpolys <- lapply(areas, function(x) which(x > minarea))
  length(unlist(bigpolys))

  # Get only the big polygons and extract them
  for(i in 1:length(bigpolys)){
    if(length(bigpolys[[i]]) >= 1 && bigpolys[[i]] >= 1){
      poly@polygons[[i]]@Polygons <- poly@polygons[[i]]@Polygons[bigpolys[[i]]]
      poly@polygons[[i]]@plotOrder <- 1:length(poly@polygons[[i]]@Polygons)
    }
  }
  return(poly)
}

provinces <- c("Ontario", "Manitoba", "Québec")

# Get ALL provincial boundaries as spatial polygons data frame
canada <- getData("GADM", country = "CAN", level = 1)
usa <- getData("GADM", country = "USA", level = 1)

# Narrow down to just the provinces of interest
ca.provinces <- canada[canada$NAME_1 %in% provinces,]

ca.provinces_small <- getSmallPolys(ca.provinces) %>% 
  gSimplify(tol = 0.01, topologyPreserve = TRUE) %>% 
  poly2df()

usa_small <- getSmallPolys(usa) %>% 
  gSimplify(tol = 0.01, topologyPreserve = TRUE) %>% 
  poly2df()


map_theme <- theme_void() +
  theme(
    panel.background = element_rect(fill = "steelblue")
  )


ggplot() +

  # Crop area
  coord_fixed(xlim = c(-95, -65),  ylim = c(41, 55), ratio = 1.6) +

  # Canada
  geom_polygon(data = ca.provinces_small,  aes(lon, lat), 
               fill = "dark grey", size = 0.3, color = "black") +
  geom_polygon(data = usa_small,  aes(lon, lat), 
               fill = "dark grey", size = 0.3, color = "black") + 
  map_theme
  • 1
    tip: use the sf package, convert those objects to sf spatial objects, and use geom_sf - it won't help with the lake problem, but there's then no need to convert to a data frame and plotting is faster. – Spacedman Jan 29 at 8:13
2

Using the lakes shapefile referenced in the other answer, you can subtract the lakes from the boundaries to get the land region for each country. Working with sf objects helps greatly too:

library(raster)
canada <- getData("GADM", country = "CAN", level = 1)
usa <- getData("GADM", country = "USA", level = 1)
provinces <- c("Ontario", "Manitoba", "Québec")
ca.provinces <- canada[canada$NAME_1 %in% provinces,]

library(sp)
library(sf)

Now get the lakes as an sf object, and convert the provinces and usa to sf objects:

lakes = st_read("./Great_Lakes.shp")

ca.provinces = st_as_sf(ca.provinces)
usa = st_as_sf(usa)

Subtract the lakes from the admin regions:

ca.land = st_difference(ca.provinces,lakes)
usa.land = st_difference(usa,lakes)

Now we use ggplot2 with geom_sf and your theme, and coord_sf to set the limits:

library(ggplot2)

ggplot() + geom_sf(data=usa.land) + geom_sf(data=ca.land) + coord_sf(xlim = c(-95, -65),  ylim = c(41, 55)) + map_theme

enter image description here

Note how the background shows through, which you don't get by simply splatting a lake polygon shape over the top with the right colour.

If things are still slow then use st_simplify to reduce the complexity of polygons. But this isn't slow on my fairly old hardware.

  • Thank-you for your clear and complete answer. – ixodid Jan 29 at 19:06
1

One solution would be download a Great Lakes shapefile (e.g. the shapes here) and add them as a last layer on your map, the same color as the background.

Edit: You can get one good Great Lakes shapefile here rather than shapes for each individual link as above.

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.