Geospatial Visualization using R

class: center, middle, inverse, title-slide

# Geospatial Visualization using R
## Part 6: Static Maps
### Bhaskar V. Karambelkar
### 2017/07/04

---

# Ways To Create Static Maps

- Base plotting

- `ggplot2` and related packages

- Ease-of-use packages

- `tmap`
  - `choroplethr`
  - `cartography`
  - `ggmap`

---

# Plotting with Base R: Sample 1

- `plot()` works on all spatial objects. `plot(..., add=TRUE)` to overlay.

- Can be cumbersome to create production-ready maps.

```r
suppressPackageStartupMessages(library(sp))
demo(meuse, ask = FALSE, echo = FALSE) 
crs.longlat = CRS("+init=epsg:4326")
meuse.longlat = spTransform(meuse, crs.longlat)
meuse.riv.longlat <- spTransform(meuse.riv, crs.longlat)
grid.lines <- gridlines(meuse.longlat)
```

```r
par(mar = c(1, 1, 1, 1))
plot(methods::as(meuse.longlat, "Spatial"), expandBB=c(0.05,0,0.1,0))
plot(grid.lines, add = TRUE, col = grey(.8))
plot(meuse.longlat, pch=1, cex = sqrt(meuse$zinc)/12, add = TRUE)
text(labels(grid.lines, side=2:3), col = grey(.7), offset=1.5)
v = c(100,200,400,800,1600)
legend("bottomright", legend = v, pch = 1, pt.cex = sqrt(v)/12,
       text.col =grey(.8), box.col=grey(0.8), title='Zinc Conc. (ppm)')
plot(meuse.riv.longlat, add = TRUE, col = grey(.9, alpha = .5))
```

---

# Plotting with Base R: Sample 1

---

# Plotting with Base R: Sample 2

```r
library(maps)
data(world.cities)
world.cities <- world.cities[world.cities$pop>1000000,]
coordinates(world.cities) <- ~long+lat 
proj4string(world.cities) <- '+init=epsg:4326'

world <- rnaturalearth::countries110
world <- world[world$name != 'Antarctica',]

grid.lines.mj <- gridlines(world,easts = seq(-180,180,by=30), norths = seq(-90,90,by=30))
grid.lines.mi <- gridlines(world,easts = seq(-165,195,by=15), norths = seq(-90,90,by=15))

par(mar = c(8, 0.1, 0.1, 0.1))
plot(methods::as(world, 'Spatial'), expandBB=c(0,0,0.05,0.05))
plot(world, add=TRUE, border=grey(0.2))
plot(grid.lines.mi, col=grey(0.95), add=T)
plot(grid.lines.mj, col=grey(0.9), add=T)
text(labels(grid.lines.mj, side=1:2, col = grey(.6), offset=0.3)
plot(world.cities, add=TRUE, col='#FF5A0088', pch=20,
     cex=world.cities$pop/2000000)
v = c(1,4,8,12)
legend("topright", legend = v, pch = 20, pt.cex = v/2,
       text.col =grey(.7), box.col=grey(0.9), col = '#FF5A0088',
       title='Pop. (Millions)', horiz =T)
```

---

# Plotting with Base R: Sample 2

---

# Plotting with Base R: Sample 2

- Actually I cheated with the last plot.

- The plot you saw as plotted in [Winkel triple](https://en.wikipedia.org/wiki/Winkel_tripel_projection) projection, while the code was using lat/long values.

- Before the `par(...)` line we need to convert all the data to our target projection.

```r
world <- spTransform(world, CRS("+proj=wintri"))
world.cities <- spTransform(world.cities, CRS("+proj=wintri"))
grid.lines.mj <- spTransform(grid.lines.mj,CRS("+proj=wintri"))
grid.lines.mi <- spTransform(grid.lines.mi,CRS("+proj=wintri"))
```

- But the lables were showing lat/long values, how come ?

```r
text(labels(grid.lines.mj, side=1:2, labelCRS = CRS("+init=epsg:4326")),
     col = grey(.6), offset=0.3)
```

- [projectionwizard.org](http://projectionwizard.org/) helps you determine sensible projection/s for your data.

---

# Plotting with Base R: Miscellaneous Notes

- `sp::gridlines()`, `sp::graiddata()`, `sp::labels()`, `sp::text()`: <br/> Creates and plot EW and NS grid lines

- `sp::bbox()` / `sp::bbexpand()`: Retrieves / expands the bounding-box of a spatial object.

- Raster data is plotted similar to vector data, but for reprojecting rasters use `raster::projectRaster()` instead of `sp::Transform()`.

- [`sp` Gallery](https://edzer.github.io/sp/) is probably the best source for advanced `sp` plotting using base R graphics.

- [`sf` Plotting](http://r-spatial.org/r/2017/01/12/newssf.html) is discussed in detail by Edzer Pebesma, the author of `sf`.

- Robin Lovelace maintains a really good [repository](https://github.com/Robinlovelace/Creating-maps-in-R) intruducing spatial visualization in R.

- Zev Ross (again!) has a very good [blog post](http://zevross.com/blog/2015/03/30/map-and-analyze-raster-data-in-r/) for working with raster data in R.

---

# Plotting with `ggplot2`: `sp` Data

```r
world <- rnaturalearth::countries110
europe <- world[world$region_un=="Europe"&world$name!='Russia',]
# plot(europe)

# Let's add a unique ID column to our data.
*europe@data$id <- row.names(europe@data)

# A bounding box for continental Europe.
europe.bbox <- SpatialPolygons(list(Polygons(list(Polygon(
  matrix(c(-25,29,45,29,45,75,-25,75,-25,29),byrow = T,ncol = 2)
)), ID = 1)), proj4string = CRS(proj4string(europe)))

# Get polygons that are only in continental Europe.
europe.clipped <-
*  rgeos::gIntersection(europe, europe.bbox, byid = TRUE, id=europe$id)
```

- `ggplot2` requires tidy data, which `sp` objects are not.

```r
europe.tidy <- broom::tidy(europe.clipped)
europe.tidy <- dplyr::left_join(europe.tidy, europe@data, by='id')
```

---

# Plotting with `ggplot2`: `sp` Data

- Ready to make a production grade map?

```r
library(ggplot2)
ggplot(europe.tidy, aes(long,lat, group=group,fill=gdp_md_est/1000)) +
  geom_polygon(alpha=0.8,color='black') +
  coord_map("azequalarea") +
  hrbrthemes::theme_ipsum_rc() +
  viridis::scale_fill_viridis(
    name='Median GDP \n(in Billions)', direction = -1, labels=scales::dollar) +
  labs(x=NULL, y=NULL, 
       title='Median GDP Estimates of\nContinental Europe & Iceland',
       caption='Source: http://www.naturalearthdata.com/')
```

---

# Plotting with `ggplot2`: `sp` Data

---

# Plotting with `ggplot2`: `sf` Data

- Similar to last plot, but with `sf`.

```r
library(sf)

world <- st_as_sf(rnaturalearth::countries110)
europe <- dplyr::filter(world, region_un=="Europe" & name!='Russia')

# A bounding box for continental Europe.
europe.bbox <- st_polygon(list(
  matrix(c(-25,29,45,29,45,75,-25,75,-25,29),byrow = T,ncol = 2)))

europe.clipped <- st_intersection(europe, st_sfc(europe.bbox, crs=st_crs(europe)))

ggplot(europe.clipped, aes(fill=gdp_md_est/1000)) +
  geom_sf(alpha=0.8,col='white') +
   coord_sf(crs="+proj=aea +lat_1=36.333333333333336 +lat_2=65.66666666666667 +lon_0=14") +
  hrbrthemes::theme_ipsum_rc() +
  viridis::scale_fill_viridis(
    name='Median GDP \n(in Billions)', direction = -1, labels=scales::dollar) +
  labs(x=NULL, y=NULL, title=NULL,
       caption='Source: http://www.naturalearthdata.com/')
```

---

# Plotting with `ggplot2`: `sf` Data

---

# Plotting with `ggplot2`:`raster` Data

```r
suppressPackageStartupMessages(library(raster))

r <- raster(system.file("external/test.grd", package="raster"))

*rasterVis::gplot(r) + geom_tile(aes(fill = value))  +
  viridis::scale_fill_viridis(direction = -1, na.value='#FFFFFF00') +
  coord_equal() + hrbrthemes::theme_ipsum()
```

---

# Plotting with `ggplot2`: Notes

- `ggplot2` is huge! `geoms`/`scales`/`legends`/`themes`/`facets`...

- There are several ggplot2 [extension](https://www.ggplot2-exts.org) packages.

- Prefer `geom_polygon` over `geom_map`.

- Be careful when tidying and subsequently joining `sp` data.

- Add the main dataset to `ggplot2()`, and others to respective `geom`'s `data` argument. You may need `inherit.aes=F` in `geom`s with their own data.

- Some good map themes `ggthemes::theme_map`, `hrbrthemes::theme_ipsum`, `ggplot2::theme_bw`.

- For `sp` objects prefer `ggalt::coord_proj` over `ggplot2::coord_map`, as the former supports proj4 strings.

- `geom_sf` support is brand new so expect a few bugs.

---

# Plotting with tmap

- `plot` mode for static maps and `view` mode for interactive maps.

- `qtm()` Quick plotting. Otherwise `tm_shape()` + `tm_<feature>`.

- Syntax similar to `ggplot2`, but returned object is not a `ggplot2` plot.

- Unlike `ggplot2` you need one `tm_shape` + one `tm_<feature>` per layer.

- Supports many different map styles and aesthetic elements like compass, legend, scale etc.

```r
usa_pop_history <- readr::read_tsv(system.file(
  'extdata','usa_pop_history.tsv', package='user2017.geodataviz'))
usa <- dplyr::left_join(albersusa::usa_sf(), usa_pop_history, by=c('name'='State'))
usa <- st_transform(usa, crs=albersusa::us_laea_proj)

years <-  c(1900,1950,1990,2000,2010,2015)

tmap::qtm(usa, fill=paste0("p.",years), title=years, 
          layout.title.color='red', layout.title.position=c('center','top'),
          layout.main.title = 'U.S.A. Population by State',
          style = 'col_blind', facets.nrow=2)
```

---

# Plotting with `tmap`

---

# Plotting with `ggmap`

- Primarily used to add base maps from Google Maps, OSM, Stamen etc. to `ggplot2` maps.

- The returned object is a `ggplot2` object which can be further worked upon.

```r
library(ggplot2)
library(ggmap)
library(magrittr)

paste0("2016-0",1:7) %>%
  purrr::map(function(month) {
    suppressMessages(readr::read_csv(
      system.file(
        sprintf("examples/data/London-Crimes/%s/%s-city-of-london-street.csv.zip",
                month,month), package='leaflet.extras')
    ))
  }) %>%
  dplyr::bind_rows() -> crimes

qmplot(Longitude, Latitude, data=crimes, geom="blank", zoom=15,
       maptype = "toner-lite", facets = ~Month) +
  stat_density_2d(aes(fill = ..level..), geom = "polygon", alpha = .3) +
  scale_fill_gradient2("Crime Heatmap", low = "white", mid = "yellow", high = "red") +
  theme(legend.position = 'hide')
```

---

# Plotting with `ggmap`

---

# Other Options for Plotting

- `cartography` package eases base R plotting quite a bit.<br/> Example code of generating a choropleth map of Internet usage by population for African countries.

```r
library(cartography)

world <- sf::st_as_sf(rnaturalearth::countries110)
internet_usage <- suppressMessages(readr::read_csv(
  'inst/extdata/africa-internet_usage-2015.csv'))

africa <- dplyr::filter(world, region_un=='Africa') %>%
  dplyr::left_join(internet_usage %>% dplyr::select(
    `Country Code`, `2015 [YR2015]`
  ) %>% dplyr::rename(iso_a3=`Country Code`, internet.usage.2015=`2015 [YR2015]`),
  by = 'iso_a3') %>%
  st_transform(crs="+proj=laea +lon_0=18.984375")

par(mar = c(0.5, 0.5, 0.5, 0.5))
plot(st_geometry(africa), border=grey(0.2), col=grey(0.9))
plot(st_centroid(africa), add=TRUE, col='black', pch=20)
*propSymbolsLayer(x=africa, var='internet.usage.2015', inches = 0.3, col = '#FF5A0088')
```

---

# Words of Wisedom

- Base R Graphics ...
  - Powerful but unintuitive.
  - `cartography` package eases certain parts.

- `tmap`
  - Extremely easy to get started and very flexible.
  - Not `ggplot2` compatible, and has its own syntax.
  - Has a `view` mode for interactive maps.

- `ggplot2`
  - Structured API and very flexible.
  - `ggmap` & `choroplethr` serve as addons providing specialized plotting functions.
  - [extensions](https://www.ggplot2-exts.org/) galore.
  - Can be made interactive using `plotly` package.

---

class: inverse middle

# Part 6: The End!

Continue to [Part 7: Interactive Maps](07-Interactive_Maps.html)

.footnote[Restart [Part 6: Static Maps](06-Static_Maps.html)]