Work code - see the slides for the final code.

cf.  - https://chiajung-yeh.github.io/Spatial-Analysis/

library(tidyverse)

## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --

## v ggplot2 3.3.5     v purrr   0.3.4
## v tibble  3.1.6     v dplyr   1.0.7
## v tidyr   1.1.4     v stringr 1.4.0
## v readr   2.1.0     v forcats 0.5.1

## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()

library(sf)

## Linking to GEOS 3.9.1, GDAL 3.2.1, PROJ 7.2.1; sf_use_s2() is TRUE

library(spData)

world is an R list that represents an simple-feature (sf) collection of 241 countries, 63 fields per country. Note:

• Bounding box of most of the globe: -180–180 logitude and -89.99–83.59 latitude
• CRS: WGS84 aka. EPSG:4326. This one, represented using a complicated string (PROJ4), has three basic components:
  • Datum: commonly WGS coordinates
  • Projection: This one is longlat which isn’t a projection - geom_sf must have a default.
  • Additional Parameters: e.g., center-of-map ?? (for Albers?)
• Geometry: MULTIPOLYGON geometry specified using:
  • Well-Known Text (WKT) format (i.e., points, linestrings, polygons, and multi-* versions of all three)
  • Coordinates: x (longitude); y (latitude) (i.e., longlat). Note that these are angular because the specify a location on the Earth’s sphere (see WGS84 reference frame.

ggplot & sf work together to plot this data easily.

world

## Simple feature collection with 177 features and 10 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -180 ymin: -89.9 xmax: 180 ymax: 83.64513
## Geodetic CRS:  WGS 84
## # A tibble: 177 x 11
##    iso_a2 name_long continent region_un subregion type  area_km2     pop lifeExp
##  * <chr>  <chr>     <chr>     <chr>     <chr>     <chr>    <dbl>   <dbl>   <dbl>
##  1 FJ     Fiji      Oceania   Oceania   Melanesia Sove~   1.93e4  8.86e5    70.0
##  2 TZ     Tanzania  Africa    Africa    Eastern ~ Sove~   9.33e5  5.22e7    64.2
##  3 EH     Western ~ Africa    Africa    Northern~ Inde~   9.63e4 NA         NA  
##  4 CA     Canada    North Am~ Americas  Northern~ Sove~   1.00e7  3.55e7    82.0
##  5 US     United S~ North Am~ Americas  Northern~ Coun~   9.51e6  3.19e8    78.8
##  6 KZ     Kazakhst~ Asia      Asia      Central ~ Sove~   2.73e6  1.73e7    71.6
##  7 UZ     Uzbekist~ Asia      Asia      Central ~ Sove~   4.61e5  3.08e7    71.0
##  8 PG     Papua Ne~ Oceania   Oceania   Melanesia Sove~   4.65e5  7.76e6    65.2
##  9 ID     Indonesia Asia      Asia      South-Ea~ Sove~   1.82e6  2.55e8    68.9
## 10 AR     Argentina South Am~ Americas  South Am~ Sove~   2.78e6  4.30e7    76.3
## # ... with 167 more rows, and 2 more variables: gdpPercap <dbl>,
## #   geom <MULTIPOLYGON [°]>

world %>%
 ggplot() +
   geom_sf()

Here’s some of the data for Afghanistan. Note that R list represents this in columns, with separate/parallel lists for: name, geometry, …

country_index <- 2  # Afghanistan

world$name[country_index]

## Warning: Unknown or uninitialised column: `name`.

## NULL

world$region_wb[country_index]

## Warning: Unknown or uninitialised column: `region_wb`.

## NULL

world$geometry[country_index]

## Warning: Unknown or uninitialised column: `geometry`.

## NULL

world$geometry[country_index] %>%
 ggplot() +
   geom_sf() +
   geom_sf_text(aes(label = world$formal_en[country_index]))

## Warning: Unknown or uninitialised column: `geometry`.

## Warning: Unknown or uninitialised column: `formal_en`.

## Warning: Unknown or uninitialised column: `formal_en`.

Because the Earth is not flat, all geospatial data needs to be projected onto a 2D surface to be visualized as a map. See the orange peel analogy.

Because the CRS for the plots above is longlat, which is unprojected, ggplot/geom_sf are using a default projection that I don’t know. The manual says: “The coordinate reference system (CRS) into which all data should be projected before plotting. If not specified, will use the CRS defined in the first sf layer of the plot.”

world %>%
  ggplot() +
    geom_sf() +
    coord_sf(crs = 4326) +  # The default, whatever it is
    ggtitle("WGS 84 - WGS84 - World Geodetic System 1984, used in GPS")

    # cowplot::theme_minimal_grid()

world %>%
  ggplot() +
    # ggspatial::layer_spatial(data = world)+
    coord_sf(crs = 3395) +
    ggtitle("WGS 84 / World Mercator")+
    cowplot::theme_minimal_grid()