Useful websites to reference for these tools:

Mapview Controls: https://r-spatial.github.io/mapview/articles/mapview_02-advanced.html

Mapview Basemaps: https://leaflet-extras.github.io/leaflet-providers/preview/

ggplot mapping: https://r-spatial.github.io/sf/articles/sf5.html

Producing kable tables: https://cran.r-project.org/web/packages/kableExtra/vignettes/awesome_table_in_html.html#Table_Styles

Setup

## Load the Packages for this

library(mapview)
library(raster)
library(RColorBrewer)
library(sf)
library(tmap)
library(viridis)
library(ncdf4)
library(leaflet) ## Basemaps

setwd("~/Desktop/University of Utah PhD /Course Work/Fall 2022 Semester/GEOG6000_Data Analysis/lab08b")

Static Maps with tmap

## tmap works similar to ggplot

##Using NY shapefile as an example

NY8 = st_read("../datafiles/NY_data/NY8_utm18.shp")
## Reading layer `NY8_utm18' from data source 
##   `/Users/davidleydet/Desktop/University of Utah PhD /Course Work/Fall 2022 Semester/GEOG6000_Data Analysis/datafiles/NY_data/NY8_utm18.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 281 features and 17 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 358241.9 ymin: 4649755 xmax: 480393.1 ymax: 4808545
## Projected CRS: WGS 84 / UTM zone 18N
## Extract the Polygons for Syracuse

syracuse = NY8[NY8$AREANAME == "Syracuse city" , ]

## Simple map using tm_borders()

syracuse.tmap = tm_shape(syracuse) + tm_borders()

syracuse.tmap

syracuse.tmap +
  tm_fill("POP8")

syracuse.tmap +
  tm_fill("POP8", 
          palette = brewer.pal(5, "GnBu"), 
          style = "quantile",
          title = "Population") +
  tm_layout(legend.title.size = 1,
          legend.text.size = 0.6,
          legend.position = c("right","bottom"),
          legend.bg.alpha = 1) 

##Adding graticules, north seeking arrow, credits, etc.

syracuse.tmap +
  tm_fill("POP8", 
          palette = brewer.pal(5, "Purples"), 
          style = "quantile",
          title = "Population") +
  tm_layout(legend.title.size = 1,
          legend.text.size = 0.6,
          legend.position = c("right","bottom"),
          legend.bg.alpha = 1) +
  tm_graticules(col = "lightgray") +
  tm_compass(position = c("left", "bottom")) +
  tm_credits("By David Leydet", position = c("right", "top"))

Point Data

## Read in the data
wna.climate = read.csv("../datafiles/WNAclimate.csv")

## Set as a simple feature --> point to coordinates (Long/Lat), and give it a coordinate reference system (CRS)
wna.climate = st_as_sf(wna.climate,
                       coords = c("LONDD", "LATDD"),
                       crs = 4326)


## Bring in the countries polygon shapefile

countries = st_read("../datafiles/ne_50m_admin_0_countries/ne_50m_admin_0_countries.shp", quiet = TRUE)


## Plot the map on a color scale using January Temperatures (Jan_Tmp)

tm_shape(wna.climate) +
  tm_symbols(col = "Jan_Tmp")

##Reverse the color palette using a "-" symbol

##Change the legend title using title.col =

tm_shape(wna.climate) +
  tm_symbols(col = "Jan_Tmp",
             alpha = 0.5,
             palette = "-RdBu",
             title.col = "January Temperature")

## Add the natural earth shape file (variable = countries) to display borders

## Ensure you set the bounding box for the WNAclimate data using the bbox = st_bbox syntax


tm_shape(countries, bbox = st_bbox(wna.climate)) +
  tm_borders(col = "gray") +
  tm_shape(wna.climate) +
  tm_symbols(col = "Jan_Tmp",
             alpha = 0.5,
             palette = "-RdBu",
             title.col = "January Temperature") +
  tm_layout(main.title = "W. North America Climate",
            legend.position = c("left", "bottom"))

## Add tm_style() to explore different map styles!

wna_map = 
  tm_shape(countries, bbox = st_bbox(wna.climate)) +
  tm_style("cobalt") +
  tm_borders(col = "gray") +
  tm_shape(wna.climate) +
  tm_symbols(col = "Jan_Tmp",
             alpha = 0.5,
             palette = "-RdBu",
             title.col = "January Temperature") +
  tm_layout(main.title = "W. North America Climate",
            legend.position = c("left", "bottom"))

wna_map

## Save the map using the tmap_save() function

tmap_save(wna_map, "wna_temp_jan.pdf")

Facetted Maps

## Useful for timeseries visualizations!

## Using built in data set metro

data(metro)
head(metro)
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## Simple feature collection with 6 features and 12 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -64.18105 ymin: -34.60508 xmax: 69.17246 ymax: 36.7525
## Geodetic CRS:  WGS 84
##            name             name_long iso_a3 pop1950 pop1960 pop1970 pop1980
## 2         Kabul                 Kabul    AFG  170784  285352  471891  977824
## 8       Algiers El Djazair  (Algiers)    DZA  516450  871636 1281127 1621442
## 13       Luanda                Luanda    AGO  138413  219427  459225  771349
## 16 Buenos Aires          Buenos Aires    ARG 5097612 6597634 8104621 9422362
## 17      Cordoba               Cordoba    ARG  429249  605309  809794 1009521
## 25      Rosario               Rosario    ARG  554483  671349  816230  953491
##     pop1990  pop2000  pop2010  pop2020  pop2030                    geometry
## 2   1549320  2401109  3722320  5721697  8279607   POINT (69.17246 34.52889)
## 8   1797068  2140577  2432023  2835218  3404575     POINT (3.04197 36.7525)
## 13  1390240  2591388  4508434  6836849 10428756   POINT (13.23432 -8.83682)
## 16 10513284 12406780 14245871 15894307 16956491 POINT (-58.40037 -34.60508)
## 17  1200168  1347561  1459268  1562509  1718192  POINT (-64.18105 -31.4135)
## 25  1083819  1152387  1298073  1453814  1606993 POINT (-60.63932 -32.94682)
## This is an sf object, with a set of columns with population size for each decade between 1950 and 2030. To plot this as a facetted map, we first need to convert this to long format with one column for area, one for population size and one for year. We’ll do this imply by concatenating the different populations sizes, and repeating the location names and year, then put it all into a temporary data frame:

## Creating a variable for just population

pop = c(metro$pop1950, metro$pop1960, metro$pop1970, metro$pop1980, metro$pop1990, metro$pop2000)

## Repeat the names six times to match the population data
name = rep(metro$name, 6)

## Repeat the years from 1950 to 2000 for each number of rows in metro
year = rep(seq(1950, 2000, by = 10), each = nrow(metro))

##Put it into a temporary dataframe
tmp = data.frame(pop, name, year)

## Merge it with the original simple feature using name

metro_long = merge(metro, tmp, by = "name")
tm_shape(metro_long) +
  tm_facets(by = "year") +
  tm_symbols("pop", 
             col = "firebrick",
             title.size = "Population",
             legend.size.is.portrait = TRUE) +
  tm_shape(countries) +
  tm_borders() +
  tm_layout(main.title = "Metropolitan Center Population")

Raster Data

## Reading in the data

air_temp = rotate(raster("../datafiles/air.mon.ltm.nc", varname = "air"))

## z-value gives you the date? 12/30? Basically January?
air_temp
## class      : RasterLayer 
## dimensions : 73, 144, 10512  (nrow, ncol, ncell)
## resolution : 2.5, 2.5  (x, y)
## extent     : -178.75, 181.25, -91.25, 91.25  (xmin, xmax, ymin, ymax)
## crs        : NA 
## source     : memory
## names      : Monthly.Long.Term.Mean.Air.Temperature.at.sigma.level.0.995 
## values     : -37.77505, 33.58953  (min, max)
## z-value    : 0000-12-30
## Change the layer name from Monthly.Long.Term.Mean.Air.Temperature.at.sigma.level.0.995 to something simple jan_tmp

names(air_temp) = "jan_tmp"

## Check CRS
crs(air_temp)
## Coordinate Reference System:
## Deprecated Proj.4 representation: NA
## Warning in wkt(x): CRS object has no comment
## Apply the CRS to it

proj4string(air_temp) = CRS("+init=epsg:4326")

## Double check to ensure the CRS was applied

crs(air_temp)
## Coordinate Reference System:
## Deprecated Proj.4 representation: +proj=longlat +datum=WGS84 +no_defs 
## WKT2 2019 representation:
## GEOGCRS["WGS 84",
##     ENSEMBLE["World Geodetic System 1984 ensemble",
##         MEMBER["World Geodetic System 1984 (Transit)",
##             ID["EPSG",1166]],
##         MEMBER["World Geodetic System 1984 (G730)",
##             ID["EPSG",1152]],
##         MEMBER["World Geodetic System 1984 (G873)",
##             ID["EPSG",1153]],
##         MEMBER["World Geodetic System 1984 (G1150)",
##             ID["EPSG",1154]],
##         MEMBER["World Geodetic System 1984 (G1674)",
##             ID["EPSG",1155]],
##         MEMBER["World Geodetic System 1984 (G1762)",
##             ID["EPSG",1156]],
##         MEMBER["World Geodetic System 1984 (G2139)",
##             ID["EPSG",1309]],
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",7030]],
##         ENSEMBLEACCURACY[2.0],
##         ID["EPSG",6326]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433],
##         ID["EPSG",8901]],
##     CS[ellipsoidal,2],
##         AXIS["longitude",east,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433,
##                 ID["EPSG",9122]]],
##         AXIS["latitude",north,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433,
##                 ID["EPSG",9122]]],
##     USAGE[
##         SCOPE["unknown"],
##         AREA["World."],
##         BBOX[-90,-180,90,180]]]

Raster Plot

tm_shape(air_temp) +
  tm_raster(col = "jan_tmp",
            n = 9,
            style = "fisher",
            palette = "-RdBu",
            title = "January Temperature",
            legend.hist = TRUE) +
  tm_shape(countries) +
  tm_borders() +
  tm_layout(legend.outside = TRUE,
            legend.outside.position = "left")

Tips and Tricks

## Use tm_tip() for helpful tips!

tmap_tip()
## Use size.max to set the value mapped to maximum symbol/bubble size. Use scale to increase (or decrease) all symbols/bubbles. 
## New since tmap 1.1 
## 
## data(metro)
## tm_shape(metro) + tm_bubbles("pop2020", size.max = 5e7, scale = 2) + tm_legend(outside = TRUE, outside.position = "bottom")
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()

Interactive Maps with mapview

## Not working correctly*******
my.pal = brewer.pal(5, "Reds")

## Initial mapview start

mapview(wna.climate,
        zcol = "Jul_Tmp", ## Color by July Temperature
        color.regions = brewer.pal(5, "Reds"), ##Not working
        color = "gray", ## Outline color
        alpha.regions = 0.2, ## Fill transparency
        alpha = 0.3) ## Outline Transparency
mapview(syracuse,
        col.regions = "darkseagreen",
        color = "black")
## Read in the Data

fn = system.file("extdata", "kiliNDVI.tif", package = "mapview" )

## Extract the first layer
kili_ndvi = raster::stack(fn)[[1]]

mapview(kili_ndvi,
        col.regions = viridis::viridis)

Basemaps

## Check out leaflet for basemap options - see https://leaflet-extras.github.io/leaflet-providers/preview/ for examples

mapview(syracuse,
        col.regions = "coral",
        color = "gray",
        map.types = "CartoDB.DarkMatter")

Layers

Syracuse.Centers = st_centroid(syracuse)

lyr1 = mapview(syracuse,
        col.regions = "coral",
        color = "gray",
        map.types = "CartoDB.DarkMatter")


lyr2 = mapview(Syracuse.Centers,
               col.regions = "dodgerblue",
               color = "gray",
               cex = 3,
               map.types = "CartoDB.DarkMatter")

lyr1 + lyr2

Aesthetic Mapping

## Use the zcol argument to map attributes

my.pal2 = brewer.pal(n=7, "GnBu")

mapview(syracuse,
        zcol = "POP8",
        col.regions = my.pal2,
        map.types = "CartoDB.Positron")
## Breweries in Franconia
fran.brew = mapview::breweries

## Set point size to number of different types of beer served using cex syntax
mapview(fran.brew,
        cex = "number.of.types",
        color = "gray",
        col.regions = "darkorange4",
        alpha.region = 0.4,
        layer.name = "Franconia Breweries - Number of Beers") ##Changes legend title

Inset Map

## Create a function to add the inset
add.inset = function(x) leaflet::addMiniMap(x@map)

## Create a variable to store the map
imap = mapview(fran.brew,
        cex = "number.of.types",
        color = "gray",
        col.regions = "goldenrod2",
        alpha.region = 0.4,
        layer.name = "Franconia Breweries - Number of Beers")

## Combine the inset function with your map variable
add.inset(imap)