Lecture 10

Draw Maps with sf and Google Maps

Byeong-Hak Choe
bchoe@geneseo.edu

SUNY Geneseo

April 14, 2025

Draw maps with sf

Shapefiles

  • Shapefiles are special types of data that include information about geography, such as points (latitude, longitude), paths (a bunch of connected latitudes and longitudes) and areas.

  • Most government agencies provide shapefiles for their jurisdictions.

  • For global mapping data, you can use:

  • Natural Earth

  • US Census Bureau

  • NY GIS Clearinghouse

Projections & Coordinate Systems

  • Projections matter a lot for maps.

  • You can convert your geographic data between different coordinate reference systems (CRS) with sf.

  • For instance:

    • coord_sf(crs = st_crs("XXXX")) (changes CRS on the fly in a plot)
    • st_transform() (permanently converts a data frame to a new CRS)

  • There are thousands of projections (epsg.io).

  • You need to specify both the catalog (e.g. "ESRI" or "EPSG") and the number (e.g. "54030" for Robinson).

Common Projections

  • ESRI:54002 = Equidistant Cylindrical (Gall-Peters)

  • EPSG:3395 = Mercator

  • ESRI:54009 = Mollweide

  • ESRI:54030 = Robinson

  • EPSG:4326 = WGS 84 (GPS system; –180 to 180)

  • EPSG:4269 = NAD 83 (common in North America)

  • ESRI:102003 = Albers for the contiguous US

  • Alternatively, you can pass projection names like "+proj=merc", "+proj=robin", etc.

Shapefiles to Download

Load and look at data

library(tidyverse)  # For ggplot, dplyr, and friends
library(sf)         # For GIS magic

# Example shapefile loading:
world_map <- read_sf("/Users/bchoe/My Drive/suny-geneseo/spring2025/maps/ne_110m_admin_0_countries/ne_110m_admin_0_countries.shp")
us_states <- read_sf("/Users/bchoe/My Drive/suny-geneseo/spring2025/maps/cb_2022_us_state_20m/cb_2022_us_state_20m.shp")
us_counties <- read_sf("/Users/bchoe/My Drive/suny-geneseo/spring2025/maps/cb_2022_us_county_5m/cb_2022_us_county_5m.shp")
us_states_hires <- read_sf("/Users/bchoe/My Drive/suny-geneseo/spring2025/maps/ne_10m_admin_1_states_provinces/ne_10m_admin_1_states_provinces.shp")
rivers_global <- read_sf("/Users/bchoe/My Drive/suny-geneseo/spring2025/maps/ne_10m_rivers_lake_centerlines/ne_10m_rivers_lake_centerlines.shp")
rivers_na <- read_sf('/Users/bchoe/My Drive/suny-geneseo/spring2025/lecture-code/maps/ne_10m_rivers_north_america/ne_10m_rivers_north_america.shp')
lakes <- read_sf("/Users/bchoe/My Drive/suny-geneseo/spring2025/maps/ne_10m_lakes/ne_10m_lakes.shp")
ny_schools <- read_sf('/Users/bchoe/My Drive/suny-geneseo/spring2025/lecture-code/maps/tl_2021_36_unsd/tl_2021_36_unsd.shp')

Basic Plotting

# Remove Antarctica for plotting
world_sans_antarctica <- world_map |> 
  filter(ISO_A3 != "ATA")

# Basic map
ggplot() + 
  geom_sf(data = world_sans_antarctica)

With color fills

ggplot() + 
  geom_sf(data = world_sans_antarctica, 
          fill = "#669438", 
          color = "#32481B", 
          size = 0.25) +
  theme_void()

Using the built-in ‘MAPCOLOR7’ from Natural Earth

ggplot() + 
  geom_sf(data = world_sans_antarctica, 
          aes(fill = 
                as.factor(MAPCOLOR7)),
          color = "#401D16", 
          size = 0.25) +
  scale_fill_viridis_d(
    option = "plasma") +
  guides(fill = "none") +
  theme_void()

World map: Different Projections

Robinson projection

ggplot() + 
  geom_sf(data = world_sans_antarctica, 
          fill = "#669438", color = "#32481B", size = 0.25) +
  coord_sf(crs = st_crs("ESRI:54030")) +  # Robinson
  theme_void()

Sinusoidal

ggplot() + 
  geom_sf(data = world_sans_antarctica, 
          fill = "#669438", color = "#32481B", size = 0.25) +
  coord_sf(crs = st_crs("ESRI:54008")) +
  theme_void()

Mercator

ggplot() + 
  geom_sf(data = world_sans_antarctica, 
          fill = "#669438", color = "#32481B", size = 0.25) +
  coord_sf(crs = st_crs("EPSG:3395")) +
  theme_void()

US map: Different Projections

US map

lower_48 <- us_states |> 
  filter(!(NAME %in% c("Alaska", "Hawaii", "Puerto Rico")))

# NAD83
ggplot() + 
  geom_sf(data = lower_48, fill = "#192DA1", color = "white", size = 0.25) +
  coord_sf(crs = st_crs("EPSG:4269")) +
  theme_void()

Albers

ggplot() + 
  geom_sf(data = lower_48, fill = "#192DA1", color = "white", size = 0.25) +
  coord_sf(crs = st_crs("ESRI:102003")) +
  theme_void()

tigris can “shift” AK/HI/PR to lower 48

library(tigris)
st_crs(us_states)
Coordinate Reference System:
  User input: NAD83 
  wkt:
GEOGCRS["NAD83",
    DATUM["North American Datum 1983",
        ELLIPSOID["GRS 1980",6378137,298.257222101,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["latitude",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["longitude",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4269]]
us_states_shifted <- us_states |> 
  shift_geometry()  ## from the `tigris` package
ggplot() +
  geom_sf(data = us_states_shifted) +
  coord_sf(crs = st_crs("ESRI:102003")) +
  theme_void()

Shift them “outside”

us_states_shifted_outside <- us_states |> 
  shift_geometry(position = "outside")

ggplot() +
  geom_sf(data = us_states_shifted_outside) +
  coord_sf(crs = st_crs("ESRI:102003")) +
  theme_void()

Individual States

only_ny <- lower_48 |> 
  filter(NAME == "New York")

ggplot() +
  geom_sf(data = only_ny, fill = "#1C4982") +
  theme_void()

Higher resolution from Natural Earth

only_ny_high <- us_states_hires |> 
  filter(iso_3166_2 == "US-NY")

ggplot() +
  geom_sf(data = only_ny_high, fill = "#1C4982") +
  theme_void() +
  coord_sf(crs = st_crs("EPSG:2239"))

Plotting Multiple Layers

Plotting Multiple Layers

ny_counties <- us_counties |> 
  filter(STATEFP == "36")

ggplot() +
  geom_sf(data = ny_counties) +
  theme_void()

Add state boundary + counties

ggplot() +
  geom_sf(data = only_ny_high, color = "#1C4982", size = 3) +
  geom_sf(data = ny_counties, fill = "#A5D46A", color = "white") +
  theme_void()

Add state boundary + counties

# Create a single unioned polygon from the counties:
only_ny_high_adjusted <- ny_counties |>
  summarize(geometry = st_union(geometry))

ggplot() +
  geom_sf(data = ny_counties, fill = "#A5D46A", color = "white") +
  geom_sf(data = only_ny_high_adjusted, fill=NA, color = "black", size = 3) +
  theme_void()

NY Counties

roc_counties <- ny_counties |> 
  filter(NAME %in% c("Monroe","Livingston","Ontario","Orleans",
                     "Genesee","Wyoming","Allegany", "Steuben"))
ggplot() +
  geom_sf(data = only_ny_high_adjusted, fill = "#1C4982") +
  geom_sf(data = roc_counties, fill = "#A5D46A", color = "white") +
  theme_void()

Intersection & Cutting Shapes

Example: Keep Just the Rivers in NY

# To ensure both data sets share the same CRS
st_crs(only_ny)
st_crs(rivers_na)

only_ny_4326 <- only_ny |> 
  st_transform(crs = st_crs("EPSG:4326"))

ny_rivers_na <- st_intersection(only_ny_4326, rivers_na)

ggplot() +
  geom_sf(data = only_ny) +
  geom_sf(data = ny_rivers_na, color = "#1C4982") +
  theme_void()

For a bigger water map

# To ensure data sets share the same CRS
st_crs(only_ny_high_adjusted)
st_crs(rivers_na)
st_crs(rivers_global)
st_crs(lakes)

only_ny_high_adjusted_crs <- st_transform(
  only_ny_high_adjusted,
  crs = st_crs(rivers_na)
)

ny_rivers_na <- st_intersection(
  only_ny_high_adjusted_crs,
  rivers_na
)
ny_rivers_global <- st_intersection(
  only_ny_high_adjusted_crs,
  rivers_global
)
lakes_crs <- lakes |> 
  st_transform(crs = st_crs(only_ny_high_adjusted)) |> 
  st_make_valid()
ny_lakes <- st_intersection(only_ny_high_adjusted, lakes_crs)

ggplot() +
  geom_sf(data = only_ny_high, size = 0.1) +
  geom_sf(data = ny_rivers_global, size = 0.3, color = "#1C4982") +
  geom_sf(data = ny_rivers_na, size = 0.15, color = "#1C4982") +
  geom_sf(data = ny_lakes, size = 0.3, fill = "#1C4982", color = NA) +
  coord_sf(crs = st_crs("EPSG:4326")) +
  theme_void()

Plotting Schools in New York

NY Schools

ggplot() +
  geom_sf(data = ny_schools)

Color by Water

ggplot() +
  geom_sf(data = ny_schools) +
  geom_sf(data = ny_schools, aes(fill = AWATER), size = 0.75, alpha = 0.5,
          show.legend = F) +
  theme_void()

Creating Your Own Geometry

ny_cities <- tribble(
  ~city,           ~lat,       ~long,
  "New York City", 40.712776,  -74.005974,
  "Albany",        42.652580,  -73.756230,
  "Buffalo",       42.886447,  -78.878369,
  "Rochester",     43.156578,  -77.608849,
  "Syracuse",      43.048122,  -76.147424
)


ny_cities_geometry <- ny_cities |> 
  st_as_sf(coords = c("long","lat"), crs = st_crs("EPSG:4326"))

ggplot() +
  geom_sf(data = only_ny_high_adjusted, fill = "#1C4982") +
  geom_sf(data = ny_cities_geometry, size = 3) +
  theme_void()

Add labels with geom_sf_label()

ggplot() +
  geom_sf(data = only_ny_high_adjusted, fill = "#1C4982") +
  geom_sf(data = ny_cities_geometry, size = 3) +
  geom_sf_label(aes(label = city), data = ny_cities_geometry, nudge_y = 0.2) +
  theme_void()

Automatic Geocoding

Example with tidygeocoder

library(tidygeocoder)
some_addresses <- tribble(
  ~name, ~address,
  "The White House", "1600 Pennsylvania Ave NW, Washington, DC",
  "State University of New York at Geneseo", "1 College Cir, Geneseo, NY 14454"
)

geocoded_addresses <- some_addresses |> 
  geocode(address, method = "census")

# Convert to sf
addresses_geometry <- geocoded_addresses |> 
  st_as_sf(coords = c("long","lat"), crs = st_crs("EPSG:4326"))

# Plot on a US map
ggplot() + 
  geom_sf(data = lower_48, fill = "#1C4982", color = "white", size = 0.25) +
  geom_sf(data = addresses_geometry, size = 5, color = "white") +
  geom_sf_label(data = addresses_geometry, aes(label = name),
                size = 4, fill = "white", nudge_y = 175000) + 
  coord_sf(crs = st_crs("ESRI:102003")) +
  theme_void()

Plotting Other Data on Maps

Getting World Bank data with WDI

# install.packages('WDI')
library(WDI)

WDIsearch('life expectancy')

indicators <- c(life_expectancy = "SP.DYN.LE00.IN")
wdi_raw <- WDI(country = "all", 
               indicators, 
               extra = TRUE, 
               start=2015, end=2015)
  • The WDI package allows users to search and download data from over 40 datasets hosted by the World Bank, including the World Development Indicators (‘WDI’), International Debt Statistics, Doing Business, Human Capital Index, and Sub-national Poverty indicators.

Example with World Bank data

# Suppose wdi_raw is loaded
wdi_clean_small <- wdi_raw |> select(life_expectancy, iso3c)

world_map_with_life_expectancy <- world_sans_antarctica |> 
  left_join(wdi_clean_small, by = c("ISO_A3" = "iso3c"))

ggplot() + 
  geom_sf(data = world_map_with_life_expectancy, aes(fill = life_expectancy), size=0.25) +
  coord_sf(crs = st_crs("ESRI:54030")) +
  scale_fill_viridis_c(option="viridis") +
  labs(fill = "Life Expectancy") +
  theme_void() +
  theme(legend.position = "bottom")

Fix for France & Norway

world_sans_antarctica_fixed <- world_sans_antarctica |> 
  mutate(ISO_A3 = case_when(
    ADMIN == "Norway" ~ "NOR",
    ADMIN == "France" ~ "FRA",
    TRUE ~ ISO_A3
  )) |> 
  left_join(wdi_clean_small, by = c("ISO_A3" = "iso3c"))

ggplot() + 
  geom_sf(data = world_sans_antarctica_fixed, aes(fill = life_expectancy), size=0.25) +
  coord_sf(crs = st_crs("ESRI:54030")) +
  scale_fill_viridis_c(option="viridis") +
  labs(fill = "Life Expectancy") +
  theme_void() +
  theme(legend.position = "bottom")

Plotting Custom Google Maps with ggmap

Getting a Google Maps API Key for ggmap

  • To use the ggmap package with Google Maps, you need an API key. Here’s how to get it:
  1. Go to Google Cloud Console
  2. Create or Select a Project
    • Click the project dropdown and choose or create a project.
  3. Enable the Maps API
    • Navigate to “APIs & Services” > “Library”
    • Search for and enable:
      • Maps JavaScript API
      • Geocoding API
      • Static Maps API (optional)
  4. Get Your API Key
    • Go to “APIs & Services” > “Credentials”
    • Click “Create credentials” > “API key”
    • Copy the generated key.

Setup: Google API Key

library(ggmap)
library(gmapsdistance)
google_api <- "YOUR_GOOGLE_API_KEY"
register_google(google_api) # for the ggmap package
set.api.key(google_api)     # for the gmapsdistance package
  • google_api: Replace “my_api” with your actual Google API key.
  • register_google(): Allows ggmap to use your Google Maps API.
  • set.api.key(): Allows gmapsdistance to connect to Google’s Distance Matrix API.

Basic Map Calls

NYC_Map <- get_map("New York City", 
                   source = "google", 
                   api_key = apiKey, 
                   zoom = 10)

ggmap(NYC_Map) +
  theme_map()
  • get_map(): Downloads a static Google map for New York City at zoom level 10.
  • ggmap(NYC_Map): Plots the map background.

Another Example

ROC_Map <- get_map("Rochester, NY", 
                   source = "google", 
                   api_key = apiKey, 
                   zoom = 9) # the larger value, the more zoom-in

ggmap(ROC_Map) +
  theme_map()
  • Same get_map() call, but for Rochester, NY.

Plotting Specific Locations

Geneseo_Map <- get_map("Geneseo, NY",
                       source = "google",
                       api_key = apiKey,
                       zoom = 14) 

locations <- geocode(
  c("South Hall, Geneseo, NY", 
    "Newton Hall, Geneseo, NY")
)

locations <- locations |>
  mutate(label = c("South Hall", "Newton Hall"))

ggmap(Geneseo_Map) + 
  geom_point(data = locations, size = 1, alpha = .5) +
  geom_text_repel(data = locations, 
                  aes(label = label),
                  box.padding = 1.75) +
  theme_map()
  • geocode(): Returns latitude and longitude of the string address.

Custom Fonts with showtext

library(showtext)
font_add_google("Roboto", "roboto")
font_add_google("Lato", "lato")
font_add_google("Poppins", "poppins")
font_add_google("Nunito", "nunito")
font_add_google("Annie Use Your Telescope", "annie")
font_add_google("Pacifico", "pacifico")
showtext_auto()  # Enables Google fonts in R plots.
  • showtext can load fonts from various sources, including system fonts, font files, and Google Fonts.
  • font_add_google(): simplifies adding Google Fonts by fetching and registering fonts directly from the Google Fonts repository using just the font family name.
  • showtext_auto() function enables automatic font rendering for all graphics.

Plot again with a custom font

ggmap(Geneseo_Map) + 
  geom_point(data = locations, size = 1, alpha = .5) +
  geom_text_repel(data = locations, 
                  aes(label = label),
                  family = "pacifico",
                  box.padding = 1.75) +
  theme_map()

Adding Routes

route_df <- route(
  from = "South Hall, Geneseo, NY",
  to   = "Newton Hall, Geneseo, NY",
  mode = "walking",
  structure = "route"  # returns points along the path 
)

ggmap(Geneseo_Map) +
  geom_path(data = route_df, 
            aes(x = lon, y = lat), 
            color = "#1C4982", size = 2, 
            lineend = "round") +  # round, butt, square
  geom_point(data = locations, 
             aes(x = lon, y = lat), 
             size = 2, color = "darkred") +
  geom_text_repel(data = locations, 
                  aes(x = lon, y = lat, label = label),
                  family = "nunito", 
                  box.padding = 1.75) +
  theme_map()
  • route(): Calculates a route between two addresses using Google Maps directions.
  • mode: Type of travel (walking, driving, bicycling, transit).
  • geom_path(): Draws the route line from route_df coordinates.
  • geom_point(): Marks the start/end points.

Different Map Types

# Available maptype options:
# "roadmap", "satellite", "terrain", "hybrid"
Geneseo_Map_Hybrid <- get_map("Geneseo, NY",
                              source = "google",
                              api_key = apiKey,
                              zoom = 14, 
                              maptype = "hybrid")

ggmap(Geneseo_Map) +
  geom_path(data = route_df, 
            aes(x = lon, y = lat), 
            color = "#1C4982", 
            size = 2, 
            lineend = "round",  # round, butt, square
            arrow = arrow() # Arrow option
            ) +  
  geom_point(data = locations, 
             aes(x = lon, y = lat), 
             size = 2, color = "darkred") +
  geom_text_repel(data = locations, 
                  aes(x = lon, y = lat, label = label),
                  family = "nunito", 
                  box.padding = 1.75) +
  theme_map()

Plotting Multiple Markers with Custom Icons

# install.packages("ggimage")
library(ggimage)

# Let's suppose we have an icon URL
icon_url <- "https://bcdanl.github.io/lec_figs/marker-icon-red.png"

ggmap(Geneseo_Map) +
  geom_image(data = locations,
             aes(x = lon, y = lat, image = icon_url),
             size = 0.05) +  # tweak size
  geom_text_repel(data = locations, 
                  aes(label = label),
                  box.padding = 1.75) +
  theme_map()