Module15_c

Load all the files

library(sf)

Linking to GEOS 3.13.1, GDAL 3.11.0, PROJ 9.6.0; sf_use_s2() is TRUE

library(terra)

terra 1.8.54

sf_use_s2(FALSE)

Spherical geometry (s2) switched off

Example 1: Salt Lake County

Make a map with vector layers (census tracts, light rail tracks, stations)

tracts <- st_read("./slc_tract/slc_tract/slc_tract_2015.shp", quiet = TRUE)
lightrail <- st_read("./LightRail_UTA/LightRail_UTA/LightRail_UTA.shp", quiet = TRUE)
stations <- st_read("./LightRailStations_UTA/LightRailStations_UTA/LightRailStations_UTA.shp", quiet = TRUE)

st_crs(tracts)$epsg

[1] 4269

st_crs(lightrail)$epsg

[1] 26912

tracts <- st_transform(tracts, st_crs(lightrail))

A simple map showing the polygon outlines

library(tmap)
tm_shape(tracts) + tm_borders()

Fill the map

tm_shape(tracts) + tm_fill("density") + tm_borders()

Fill the map with color palette

tm_shape(tracts) + tm_fill("density", palette = "Greens", style = "quantile") + tm_borders("lightgray")

── 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>)'
[cols4all] color palettes: use palettes from the R package cols4all. Run
`cols4all::c4a_gui()` to explore them. The old palette name "Greens" is named
"brewer.greens"
Multiple palettes called "greens" found: "brewer.greens", "matplotlib.greens". The first one, "brewer.greens", is returned.

Reversed color scale

tm_shape(tracts) + tm_fill("density", palette = "-magma", style = "cont") + tm_borders("lightgray")

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

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

Overlay the light rail track on the map

tm_shape(tracts) + tm_fill("density", palette = "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>)'
[cols4all] color palettes: use palettes from the R package cols4all. Run
`cols4all::c4a_gui()` to explore them. The old palette name "Greens" is named
"brewer.greens"
Multiple palettes called "greens" found: "brewer.greens", "matplotlib.greens". The first one, "brewer.greens", is returned.

Color the tracks

tm_shape(tracts) + tm_fill("density", palette = "Greens", style = "quantile") + tm_borders("lightgray") + tm_shape(lightrail) + tm_lines(lwd = 4, col = "ROUTE")

── 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>)'
[cols4all] color palettes: use palettes from the R package cols4all. Run
`cols4all::c4a_gui()` to explore them. The old palette name "Greens" is named
"brewer.greens"
Multiple palettes called "greens" found: "brewer.greens", "matplotlib.greens". The first one, "brewer.greens", is returned.

Removed the polygon

tm_shape(tracts) + tm_borders("lightgray") + tm_shape(lightrail) + tm_lines(lwd = 4, col = "ROUTE")

Add the stations

tm_shape(tracts) + tm_fill("density", palette = "Greens", style = "quantile") + tm_borders("lightgray") + tm_shape(lightrail) + tm_lines(lwd = 2) + tm_shape(stations) + tm_dots(size = 0.25, shape = 23)

── 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>)'
[cols4all] color palettes: use palettes from the R package cols4all. Run
`cols4all::c4a_gui()` to explore them. The old palette name "Greens" is named
"brewer.greens"
Multiple palettes called "greens" found: "brewer.greens", "matplotlib.greens". The first one, "brewer.greens", is returned.

Add details

tm_shape(tracts) + 
  tm_graticules(col = "lightgray") + 
  tm_fill("density", title = "Popn density", palette = "Greens", style = "quantile") +
  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_scale_bar(position = c("right", "top")) +
  tm_layout(main.title = "Salt Lake County Light Rail", 
            legend.outside = TRUE,
            legend.outside.position = c("left"))

── 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>)'
[v3->v4] `tm_fill()`: migrate the argument(s) related to the legend of the
visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'
! `tm_scale_bar()` is deprecated. Please use `tm_scalebar()` instead.
[v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(main.title = )`
[cols4all] color palettes: use palettes from the R package cols4all. Run
`cols4all::c4a_gui()` to explore them. The old palette name "Greens" is named
"brewer.greens"
Multiple palettes called "greens" found: "brewer.greens", "matplotlib.greens". The first one, "brewer.greens", is returned.

Clipped to smaller regions

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 ───────────────────────────────────────────────────────

Make interactive map

tmap_mode("view")

ℹ tmap mode set to "view".

tm_shape(tracts) + 
  tm_fill("density", title = "Popn density", palette = "Greens", style = "quantile", id = "TRACTCE") +
  tm_borders("lightgray") +
  tm_shape(lightrail) +
  tm_lines(lwd = 2) +
  tm_shape(stations) +
  tm_dots(size = 0.25, shape = 23)

── 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>)'[v3->v4] `tm_fill()`: migrate the argument(s) related to the legend of the
visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'[cols4all] color palettes: use palettes from the R package cols4all. Run
`cols4all::c4a_gui()` to explore them. The old palette name "Greens" is named
"brewer.greens"Multiple palettes called "greens" found: "brewer.greens", "matplotlib.greens". The first one, "brewer.greens", is returned.

tmap_mode("plot")

ℹ tmap mode set to "plot".

Example 2: Western North America climate

Convert the data frame

wna_climate <- read.csv("WNAclimate.csv")
wna_climate <- st_as_sf(wna_climate, coords = c("LONDD", "LATDD"), crs = 4326)

Make a thematic map

tm_shape(wna_climate) + tm_symbols(col = "Jul_Tmp")

Download country outlines and river centerlines

library(rnaturalearth)
countries50 <- ne_download(scale = 50, type = "admin_0_countries")

Reading layer `ne_50m_admin_0_countries' from data source 
  `C:\Users\kuma3\AppData\Local\Temp\RtmpMhFXC1\ne_50m_admin_0_countries.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 242 features and 168 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -180 ymin: -89.99893 xmax: 180 ymax: 83.59961
Geodetic CRS:  WGS 84

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

Reading layer `ne_50m_rivers_lake_centerlines' from data source 
  `C:\Users\kuma3\AppData\Local\Temp\RtmpMhFXC1\ne_50m_rivers_lake_centerlines.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 478 features and 36 fields (with 1 geometry empty)
Geometry type: MULTILINESTRING
Dimension:     XY
Bounding box:  xmin: -165.2439 ymin: -50.24014 xmax: 176.3258 ymax: 73.3349
Geodetic CRS:  WGS 84

Put this all together as a series of layers

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 = "-RdBu", alpha = 0.75, style = "cont", title.col = "degC")

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

[v3->v4] `symbols()`: instead of `style = "cont"`, use fill.scale =
`tm_scale_continuous()`.
ℹ Migrate the argument(s) 'palette' (rename to 'values') to
  'tm_scale_continuous(<HERE>)'
[v3->v4] `symbols()`: use 'fill' for the fill color of polygons/symbols
(instead of 'col'), and 'col' for the outlines (instead of 'border.col').
[v3->v4] `symbols()`: use `fill_alpha` instead of `alpha`.
[v3->v4] `symbols()`: migrate the argument(s) related to the legend of the
visual variable `fill` namely 'title.col' (rename to 'title') to 'fill.legend =
tm_legend(<HERE>)'
Multiple palettes called "rd_bu" found: "brewer.rd_bu", "matplotlib.rd_bu". The first one, "brewer.rd_bu", is returned.

[cols4all] color palettes: use palettes from the R package cols4all. Run
`cols4all::c4a_gui()` to explore them. The old palette name "-RdBu" is named
"rd_bu" (in long format "brewer.rd_bu")

Make a two panel plot maps

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 = "-RdBu", alpha = 0.75, style = "cont", title.col = "degC")

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

[v3->v4] `symbols()`: instead of `style = "cont"`, use fill.scale =
`tm_scale_continuous()`.
ℹ Migrate the argument(s) 'palette' (rename to 'values') to
  'tm_scale_continuous(<HERE>)'
[v3->v4] `symbols()`: use `fill_alpha` instead of `alpha`.
[v3->v4] `symbols()`: migrate the argument(s) related to the legend of the
visual variable `fill` namely 'title.col' (rename to 'title') to 'fill.legend =
tm_legend(<HERE>)'

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 = "BuPu", alpha = 0.75, style = "cont", title.col = "mm/yr")

[v3->v4] `symbols()`: use `fill_alpha` instead of `alpha`.
[v3->v4] `symbols()`: migrate the argument(s) related to the legend of the
visual variable `fill` namely 'title.col' (rename to 'title') to 'fill.legend =
tm_legend(<HERE>)'

tmap_arrange(m1, m2)

Multiple palettes called "rd_bu" found: "brewer.rd_bu", "matplotlib.rd_bu". The first one, "brewer.rd_bu", is returned.
[cols4all] color palettes: use palettes from the R package cols4all. Run
`cols4all::c4a_gui()` to explore them. The old palette name "-RdBu" is named
"rd_bu" (in long format "brewer.rd_bu")[cols4all] color palettes: use palettes from the R package cols4all. Run
`cols4all::c4a_gui()` to explore them. The old palette name "BuPu" is named
"brewer.bu_pu"Multiple palettes called "bu_pu" found: "brewer.bu_pu", "matplotlib.bu_pu". The first one, "brewer.bu_pu", is returned.

Example 3: Raster images

Load the data

filenames <- paste0('./rs/rs/LC08_044034_20170614_B', 1:11, ".tif")
landsat <- rast(filenames)
ca_places <- st_read("./ca_places/ca_places/ca_places.shp", quiet = TRUE)

Remake the NDVI layer

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()

[scale] tm_raster:() the data variable assigned to 'col' contains positive and negative values, so midpoint is set to 0. Set 'midpoint = NA' in 'fill.scale = tm_scale_intervals(<HERE>)' to use all visual values (e.g. colors)

Update this map

tm_shape(clamp(ndvi, lower = 0)) + tm_raster(palette = "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 = )`
[cols4all] color palettes: use palettes from the R package cols4all. Run
`cols4all::c4a_gui()` to explore them. The old palette name "Greens" is named
"brewer.greens"

Make the 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`
[v3->v4] `tm_rgb()`: instead of using r = 4, g = 3 and b = 2 , please use col =
tm_vars(c(4, 3, 2), multivariate = TRUE)

Make it a lot brighter

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`
[v3->v4] `tm_rgb()`: instead of using r = 4, g = 3 and b = 2 , please use col =
tm_vars(c(4, 3, 2), multivariate = TRUE)

Leaflet

Load library

library(leaflet)
library(htmltools)

Drop a marker on the University of Utah

m <- leaflet() %>%
  addTiles() %>%
  addMarkers(lng=-111.8421, lat=40.7649, popup="The University of Utah")
m

ESRI map

m <- leaflet() %>%
  addProviderTiles("Esri.WorldStreetMap") %>%  
  addMarkers(lng=-111.8421, lat=40.7649, popup="The University of Utah")
m

ESRI imagery

m <- leaflet() %>%
  addProviderTiles("Esri.WorldImagery") %>%  # Add default OpenStreetMap map tiles
  addMarkers(lng=-111.8421, lat=40.7649, popup="The University of Utah")
m

Add the census tracts and light rail

Reproject to WGS84 long-lat (EPSG code 4326)

tracts_ll <- st_transform(tracts, crs = 4326)
stations_ll <- st_transform(stations, crs = 4326)
lightrail_ll <- st_transform(lightrail, crs = 4326)

Add these to a dark basemap

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 some extra details

Get the list of unique station names

unique(stations_ll$LINENAME)

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

Convert the station names to an integer code

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

Make a popup window for each station

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
  
)

Rebuild the map

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
  )