Module 15 C, Better Maps

Example 1: Salt Lake County

# Load vector layers from the working directory
tracts = st_read("slc_tract_2015.shp", quiet = TRUE)
lightrail = st_read("LightRail_UTA.shp", quiet = TRUE)
stations = st_read("LightRailStations_UTA.shp", quiet = TRUE)

# Check coordinate reference systems
st_crs(tracts)$epsg

[1] 4269

st_crs(lightrail)$epsg

[1] 26912

# Reproject the tracts to match the light rail data
tracts = st_transform(tracts, st_crs(lightrail))

# Set static tmap mode
tmap_mode("plot")

ℹ tmap modes "plot" - "view"
ℹ toggle with `tmap::ttm()`

# Basic polygon outline map
tm_shape(tracts) +
  tm_borders()

# Fill polygons by density
tm_shape(tracts) +
  tm_polygons(
    fill = "density",
    fill.scale = tm_scale_continuous(values = "brewer.greens"),
    fill.legend = tm_legend(title = "Popn density")
  )

# Use a color palette with quantile breaks
tm_shape(tracts) +
  tm_polygons(
    fill = "density",
    fill.scale = tm_scale_intervals(
      values = "brewer.greens",
      style = "quantile"
    ),
    fill.legend = tm_legend(title = "Popn density")
  ) +
  tm_borders("lightgray")

# Different palette
# run cols4all::c4a_palettes() to see valid colors
tm_shape(tracts) +
  tm_polygons(
    fill = "density",
    fill.scale = tm_scale_continuous(values = "bivario.plum_mint"),
    fill.legend = tm_legend(title = "Popn density")
  ) +
  tm_borders("lightgray")

#Add light rail
tm_shape(tracts) + 
  tm_fill("density", palette = "brewer.greens", style = "quantile") +
  tm_borders("lightgray") +
  tm_shape(lightrail) +
  tm_lines(lwd = 2, lty = 'dashed', col = "darkorange")

── tmap v3 code detected ───────────────────────────────────────────────────────
[v3->v4] `tm_fill()`: instead of `style = "quantile"`, use fill.scale =
`tm_scale_intervals()`.
ℹ Migrate the argument(s) 'style', 'palette' (rename to 'values') to
  'tm_scale_intervals(<HERE>)'

# Use the route variable to theme the light rail lines
tm_shape(tracts) +
  tm_borders("lightgray") +
  tm_shape(lightrail) +
  tm_lines(lwd = 4, col = "ROUTE")

# Add the light rail tracks and stations
tm_shape(tracts) +
  tm_polygons(
    fill = "density",
    fill.scale = tm_scale_intervals(values = "brewer.greens", style = "quantile"),
    fill.legend = tm_legend(title = "Popn density")
  ) +
  tm_borders("lightgray") +
  tm_shape(lightrail) +
  tm_lines(lwd = 2) +
  tm_shape(stations) +
  tm_dots(size = 0.25, shape = 23)

# Add map elements
tm_shape(tracts) +
  tm_graticules(col = "lightgray") +
  tm_polygons(
    fill = "density",
    fill.scale = tm_scale_intervals(values = "brewer.greens", style = "quantile"),
    fill.legend = tm_legend(title = "Popn density")
  ) +
  tm_borders("lightgray") +
  tm_shape(lightrail) +
  tm_lines(lwd = 2) +
  tm_shape(stations) +
  tm_dots(size = 0.25, shape = 23) +
  tm_compass(position = c("left", "bottom")) +
  tm_scalebar(position = c("right", "top")) +
  tm_title("Salt Lake County Light Rail")

# Clip the map to downtown Salt Lake using the bbox argument
tracts_sub = tracts %>%
  dplyr::filter(TRACTCE %in% c(102500, 102600, 114000))

tm_shape(tracts, bbox = st_bbox(tracts_sub)) +
  tm_borders("lightgray") +
  tm_shape(lightrail) +
  tm_lines(lwd = 2) +
  tm_shape(stations) +
  tm_dots(size = 0.25, shape = 23) +
  tm_text("STATIONNAM", ymod = -1, bg.color = "white", size = 0.8)

── tmap v3 code detected ───────────────────────────────────────────────────────

# Set interactive tmap mode
tmap_mode("view")

ℹ tmap modes "plot" - "view"

tm_shape(tracts) +
  tm_polygons(
    fill = "density",
    fill.scale = tm_scale_intervals(values = "brewer.greens", style = "quantile"),
    fill.legend = tm_legend(title = "Popn density")
  ) +
  tm_borders("lightgray") +
  tm_shape(lightrail) +
  tm_lines(lwd = 2) +
  tm_shape(stations) +
  tm_dots(size = 0.25, shape = 23)

# Reset to static mode for the remaining examples
tmap_mode("plot")

ℹ tmap modes "plot" - "view"

Example 2: Western North America Climate

# Convert the climate table to an sf object using longitude and latitude columns
wna_climate = read.csv("WNAclimate.csv")

wna_climate = st_as_sf(
  wna_climate,
  coords = c("LONDD", "LATDD"),
  crs = 4326
)

# Make a thematic map of July temperature
tm_shape(wna_climate) +
  tm_symbols(col = "Jul_Tmp")

# Natural Earth
countries50 = ne_download(
  scale = 50,
  type = "admin_0_countries"
)

rivers50 = ne_download(
  scale = 50,
  type = "rivers_lake_centerlines",
  category = "physical"
)

# Layer country outlines, rivers, and climate points
tm_shape(countries50, bbox = st_bbox(wna_climate)) +
  tm_borders() +
  tm_shape(rivers50) +
  tm_lines("lightblue", lwd = 2) +
  tm_shape(wna_climate) +
  tm_symbols(
    col = "Jul_Tmp",
    palette = "-brewer.rd_bu",
    fill_alpha = 0.75,
    style = "cont",
    title.col = "degC"
  )

# Two-panel map
m1 = tm_shape(countries50, bbox = st_bbox(wna_climate)) +
  tm_borders() +
  tm_shape(rivers50) +
  tm_lines("lightblue", lwd = 2) +
  tm_shape(wna_climate) +
  tm_symbols(
    col = "Jul_Tmp",
    palette = "-brewer.rd_bu",
    fill_alpha = 0.75,
    style = "cont",
    title.col = "degC"
  )

m2 = tm_shape(countries50, bbox = st_bbox(wna_climate)) +
  tm_borders() +
  tm_shape(rivers50) +
  tm_lines("lightblue", lwd = 2) +
  tm_shape(wna_climate) +
  tm_symbols(
    col = "annp",
    palette = "brewer.bu_pu",
    fill_alpha = 0.75,
    style = "cont",
    title.col = "mm/yr"
  )

tmap_arrange(m1, m2)

Example 3: Raster Images

filenames = paste0("LC08_044034_20170614_B", 1:11, ".tif")
landsat = rast(filenames)

ca_places = st_read("ca_places.shp", quiet = TRUE)

# Calculate NDVI
b4 = landsat[[4]]
b5 = landsat[[5]]

ndvi = (b5 - b4) / (b5 + b4)

plot(
  ndvi,
  col = rev(terrain.colors(10)),
  main = "NDVI"
)

tm_shape(ndvi) +
  tm_raster()

[`tm_scale_intervals()`] Variable(s) "col" contains positive and negative
values, so midpoint is set to 0. Set midpoint = NA to show the full range of
visual values.
This message is displayed once per session.

tm_shape(clamp(ndvi, lower = 0)) +
  tm_raster(palette = "brewer.greens", style = "cont", title = "NDVI") +
  tm_shape(ca_places) +
  tm_borders(lwd = 2) +
  tm_layout(main.title = "Landsat 8 (2017/06/14)", 
            legend.position = c("left", "top"), 
            legend.bg.color = "white", legend.bg.alpha = 0.7)

── tmap v3 code detected ───────────────────────────────────────────────────────
[v3->v4] `tm_raster()`: instead of `style = "cont"`, use col.scale =
`tm_scale_continuous()`.
ℹ Migrate the argument(s) 'palette' (rename to 'values') to
  'tm_scale_continuous(<HERE>)'[v3->v4] `tm_raster()`: migrate the argument(s) related to the legend of the
visual variable `col` namely 'title' to 'col.legend = tm_legend(<HERE>)'[v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(main.title = )`

# A true color composite 
tm_shape(landsat) +
  tm_rgb(r = 4, g = 3, b = 2, max.value = 1)

── tmap v3 code detected ───────────────────────────────────────────────────────
[v3->v4] `tm_rgb()`: use `col.scale = tm_scale_rgb(max_color_value = 1)`
instead of `max.value = 1`

Warning in blend_grobs(grb, a$blend, dst = dst): Compositing operator '2' is
not supported by the current graphics device. Falling back to no blending. Try
png(type = "cairo") or svg().

# Stretch the Landsat values and replot the true color composite
landsat_stretch = stretch(
  landsat,
  minv = 0,
  maxv = 1,
  minq = 0.01,
  maxq = 0.99
)

tm_shape(landsat_stretch) +
  tm_rgb(r = 4, g = 3, b = 2, max.value = 1)

── tmap v3 code detected ───────────────────────────────────────────────────────
[v3->v4] `tm_rgb()`: use `col.scale = tm_scale_rgb(max_color_value = 1)`
instead of `max.value = 1`

Warning in blend_grobs(grb, a$blend, dst = dst): Compositing operator '2' is
not supported by the current graphics device. Falling back to no blending. Try
png(type = "cairo") or svg().

Leaflet

# A simple Leaflet example with default tiles and a marker
m = leaflet() %>%
  addTiles() %>%
  addMarkers(lng = -111.8421, lat = 40.7649, popup = "The University of Utah")
m

# A simple Leaflet example with an ESRI basemap
m = leaflet() %>%
  addProviderTiles("Esri.WorldStreetMap") %>%
  addMarkers(lng = -111.8421, lat = 40.7649, popup = "The University of Utah")
m

# A simple Leaflet example with imagery
m = leaflet() %>%
  addProviderTiles("Esri.WorldImagery") %>%
  addMarkers(lng = -111.8421, lat = 40.7649, popup = "The University of Utah")
m

# Salt Lake County Leaflet map
tracts_ll = st_transform(tracts, crs = 4326)
stations_ll = st_transform(stations, crs = 4326)
lightrail_ll = st_transform(lightrail, crs = 4326)

# Add the tracts, light rail tracks, and stations
leaflet() %>%
  addProviderTiles("CartoDB.DarkMatter") %>%
  addPolygons(
    data = tracts_ll,
    color = "#E2E2E2",
    opacity = 1,
    weight = 1,
    fillOpacity = 0.2
  ) %>%
  addPolylines(
    data = lightrail_ll
  ) %>%
  addMarkers(
    data = stations_ll,
    label = ~htmlEscape(STATIONNAM)
  )

# Add station colors and popup text
unique(stations_ll$LINENAME)

[1] "University"            "N/S TRAX"              "Med Center"           
[4] "Intermodal Hub"        "Airport"               "Draper"               
[7] "Mid-Jordan"            "West Valley"           "Sugar House Streetcar"

line_pal = RColorBrewer::brewer.pal(9, "Set1")
line_no = as.numeric(as.factor(stations_ll$LINENAME))
stations_ll$LINECOL = line_pal[line_no]

station_popup = paste0(
  "<b>Station: </b>", stations_ll$STATIONNAM,
  "<br>",
  "<b>Line Name: </b>", stations_ll$LINENAME,
  "<br>",
  "<b>Park n Ride: </b>", stations_ll$PARKNRIDE,
  "<br>",
  "<b>Address: </b>", stations_ll$ADDRESS
)

# Final Leaflet map with colored station markers and popups
leaflet() %>%
  addProviderTiles("Esri.WorldStreetMap") %>%
  addPolygons(
    data = tracts_ll,
    color = "#E2E2E2",
    opacity = 1,
    weight = 1,
    fillOpacity = 0.2
  ) %>%
  addPolylines(
    data = lightrail_ll
  ) %>%
  addCircleMarkers(
    data = stations_ll,
    color = ~LINECOL,
    label = ~htmlEscape(STATIONNAM),
    popup = station_popup
  )