+ - 0:00:00
Notes for current slide
Notes for next slide

Spatial Visualization

Thiyanga S. Talagala

1 / 73

R packages for spatial data analysis

  • rgdal

  • sp

  • rgeos

  • raster

  • sf

  • tmap

  • leaflet

  • ggmap

  • maptools

  • gstat

  • spatstat

  • stars

  • geosphere

  • RgoogleMaps

  • rasterVis

2 / 73

Geospatial vector data structures

Point: Individual longitude and latitude of locations

Eg: Locations, Buildings

Lines: Two or more vertices or points that are connected

Eg: Roads, rivers

Polygons: Three or more vertices and closed

Eg: Area of a country, state, district

3 / 73

What do we want to plot? Globe - 3D

4 / 73

What do we want to plot? Globe - 3D

Where do we plot? Computer screen or paper - 2D

4 / 73

Challenge

How to transform this three dimensional angular system to a two dimensional cartesian system?

5 / 73

Challenge

How to transform this three dimensional angular system to a two dimensional cartesian system?

Solution : Spatial projection

5 / 73

Types of map projections

  • Cylindrical projections

  • Conical projections

  • Planar projections

Visit here to see the visual illustration.

6 / 73

Coordinate Reference System (CRS)

Different types of CRSs:

  • Mercator

  • UTM

  • Robinson

  • Lambert

  • Sinusoidal

  • Alber

7 / 73

8 / 73

9 / 73

10 / 73

11 / 73

Why some areas are stretched and some areas are compressed?

Projecting data from a 3D round surface onto a 2D flat surface results in visual modifications to the data.

12 / 73

Why some areas are stretched and some areas are compressed?

Projecting data from a 3D round surface onto a 2D flat surface results in visual modifications to the data.

Why multiple CRS?

To optimize to best represent the:

shape

distance

area

direction

12 / 73

Shapefiles

  • Spatial data in vector format are often stored in a shapefile format

  • The structure of points, lines, and polygons are different, each individual shapefile can only contain one vector type (all points, all lines or all polygons).

  • You will not find a mixture of point, line and polygon objects in a single shapefile.

  • Shapefiles contain geographic vector data that is expressed in a particular CRS

13 / 73

  • A shapefile is made up of three required files with the same prefix ('spatial-data' in this case) but different extensions:

    spatial-data.shp: main file that stores records of each shape geometries

    spatial-data.shx: index of how the geometries in the main file relate to one-another

    spatial-data.dbf: attributes of each record

14 / 73

CRS

  • geodetic datum (e.g, WGS84)

  • the type of map projection (e.g., Mercator)

  • the parameters of the projection (e.g., location of the origin)
15 / 73

Two types of maps

  1. Simple feature map

  2. Polygon maps

16 / 73

Simple feature map

17 / 73

Example

# install.packages("devtools")
devtools::install_github("thiyangt/ceylon")
library(ceylon)
library(tidyverse)
library(sp)
library(viridis)
data(sf_sl_0)
class(sf_sl_0)
## [1] "sf"         "tbl_df"     "tbl"        "data.frame"
# work with spatial data; sp package will load with rgdal.
library(rgdal)
## Please note that rgdal will be retired during October 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
## See https://r-spatial.org/r/2023/05/15/evolution4.html and https://github.com/r-spatial/evolution
## rgdal: version: 1.6-7, (SVN revision 1203)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.4.2, released 2022/03/08
## Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/4.2/Resources/library/rgdal/gdal
## GDAL binary built with GEOS: FALSE 
## Loaded PROJ runtime: Rel. 8.2.1, January 1st, 2022, [PJ_VERSION: 821]
## Path to PROJ shared files: /Library/Frameworks/R.framework/Versions/4.2/Resources/library/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.6-1
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.
# for metadata/attributes- vectors or rasters
library(raster)
## 
## Attaching package: 'raster'
## The following object is masked from 'package:dplyr':
## 
##     select
crs(sf_sl_0)
## [1] "PROJCRS[\"SLD99 / Sri Lanka Grid 1999\",\n    BASEGEOGCRS[\"SLD99\",\n        DATUM[\"Sri Lanka Datum 1999\",\n            ELLIPSOID[\"Everest 1830 (1937 Adjustment)\",6377276.345,300.8017,\n                LENGTHUNIT[\"metre\",1]]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        ID[\"EPSG\",5233]],\n    CONVERSION[\"Sri Lanka Grid 1999\",\n        METHOD[\"Transverse Mercator\",\n            ID[\"EPSG\",9807]],\n        PARAMETER[\"Latitude of natural origin\",7.00047152777778,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8801]],\n        PARAMETER[\"Longitude of natural origin\",80.7717130833333,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8802]],\n        PARAMETER[\"Scale factor at natural origin\",0.9999238418,\n            SCALEUNIT[\"unity\",1],\n            ID[\"EPSG\",8805]],\n        PARAMETER[\"False easting\",500000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8806]],\n        PARAMETER[\"False northing\",500000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8807]]],\n    CS[Cartesian,2],\n        AXIS[\"(E)\",east,\n            ORDER[1],\n            LENGTHUNIT[\"metre\",1]],\n        AXIS[\"(N)\",north,\n            ORDER[2],\n            LENGTHUNIT[\"metre\",1]],\n    USAGE[\n        SCOPE[\"unknown\"],\n        AREA[\"Sri Lanka - onshore\"],\n        BBOX[5.86,79.64,9.88,81.95]],\n    ID[\"EPSG\",5235]]"
18 / 73
sf_sl_0
## Simple feature collection with 1 feature and 1 field
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 362203.3 ymin: 380301.9 xmax: 621918.1 ymax: 813560.9
## Projected CRS: SLD99 / Sri Lanka Grid 1999
## # A tibble: 1 × 2
##                                                                 geometry COUNTRY
##                                                       <MULTIPOLYGON [m]> <chr>  
## 1 (((481925.5 381353.7, 481922.9 381350.3, 481919 381348.2, 481914.6 38… SRI LA…
19 / 73

2. Visualizing using shapefiles: map of Sri Lanka

ggplot(sf_sl_0) +
  geom_sf()

20 / 73
ggplot(sf_sl_0) +
  geom_sf(fill='beige') +
theme_minimal()

21 / 73
library(knitr)
sf_sl_0 %>% kable()
geometry COUNTRY
MULTIPOLYGON (((481925.5 38... SRI LANKA
22 / 73

Your turn: Draw district and province maps

data:

province
district
23 / 73

Province

ggplot(data=province) +
  geom_sf()

District

ggplot(data=district) +
  geom_sf()

24 / 73
province
## Simple feature collection with 9 features and 3 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 362203.3 ymin: 380301.9 xmax: 621918.1 ymax: 813560.9
## Projected CRS: SLD99 / Sri Lanka Grid 1999
## # A tibble: 9 × 4
##                                              geometry PROVINCE Status population
## *                                  <MULTIPOLYGON [m]> <chr>    <chr>       <dbl>
## 1 (((498211.2 611042.6, 498401.7 610897.1, 498415.9 … CENTRAL  Provi…    2781000
## 2 (((609877 559315.9, 609857.5 559313, 609845.2 5593… EASTERN  Provi…    1746000
## 3 (((501928.5 712099.6, 501970.4 712098.8, 502020 71… NORTH C… Provi…    1386000
## 4 (((393087.5 629959, 393098.1 629956.7, 393103.2 62… NORTH W… Provi…    2563000
## 5 (((405858.7 700083.7, 405855.3 700082, 405851.5 70… NORTHERN Provi…    1152000
## 6 (((432117.3 521875.5, 432083.5 521852.9, 432041.7 … SABARAG… Provi…    2070000
## 7 (((481925.5 381353.7, 481922.9 381350.3, 481919 38… SOUTHERN Provi…    2669000
## 8 (((522825 568220.9, 522872.6 568217.7, 522926.1 56… UVA      Provi…    1387000
## 9 (((411888.2 438189.4, 411886.4 438182.3, 411881.7 … WESTERN  Provi…    6165000
25 / 73
district
## Simple feature collection with 26 features and 3 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 362203.3 ymin: 380301.9 xmax: 621918.1 ymax: 813560.9
## Projected CRS: SLD99 / Sri Lanka Grid 1999
## # A tibble: 26 × 4
##                                              geometry DISTRICT Status population
##  *                                 <MULTIPOLYGON [m]> <chr>    <chr>       <dbl>
##  1 (((392952.6 594064.2, 392985.8 593971.5, 392986.2… [UNKNOW… <NA>           NA
##  2 (((551311.8 581603.6, 551324.5 581599.7, 551334.5… AMPARA   Distr…     736000
##  3 (((501928.5 712099.6, 501970.4 712098.8, 502020 7… ANURADH… Distr…     943000
##  4 (((522825 568220.9, 522872.6 568217.7, 522926.1 5… BADULLA  Distr…     886000
##  5 (((609877 559315.9, 609857.5 559313, 609845.2 559… BATTICA… Distr…     579000
##  6 (((400911.5 497888.4, 401032.5 497869.3, 401035.2… COLOMBO  Distr…    2455000
##  7 (((425038.8 404072.5, 425030.7 404072.5, 425022 4… GALLE    Distr…    1135000
##  8 (((430775.7 534473.9, 430789.2 534441.8, 430790.2… GAMPAHA  Distr…    2423000
##  9 (((592923.3 452679, 592949.2 452679, 592989.8 452… HAMBANT… Distr…     668000
## 10 (((362968.6 763698.5, 362883.6 763686.9, 362787.9… JAFFNA   Distr…     621000
## # ℹ 16 more rows
26 / 73

Guess the output.

ggplot() +
  geom_sf(data=province, lwd=2, col="black") + 
  geom_sf(data=district, linetype=21, col="red")
27 / 73

Your turn:

Colour provinces according to the population.

28 / 73

Annotations

ggplot(sf_sl_0) +
  geom_sf(fill='beige') +
ggspatial::annotation_north_arrow(location = "br")+
  ggspatial::annotation_scale(location = "bl")

29 / 73

More work

library(spData)
## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`
## 
## Attaching package: 'spData'
## The following object is masked _by_ '.GlobalEnv':
## 
##     world
30 / 73
data("us_states") # from spData package
us_states
## Simple feature collection with 49 features and 6 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -124.7042 ymin: 24.55868 xmax: -66.9824 ymax: 49.38436
## Geodetic CRS:  NAD83
## First 10 features:
##    GEOID        NAME   REGION             AREA total_pop_10 total_pop_15
## 1     01     Alabama    South 133709.27 [km^2]      4712651      4830620
## 2     04     Arizona     West 295281.25 [km^2]      6246816      6641928
## 3     08    Colorado     West 269573.06 [km^2]      4887061      5278906
## 4     09 Connecticut Norteast  12976.59 [km^2]      3545837      3593222
## 5     12     Florida    South 151052.01 [km^2]     18511620     19645772
## 6     13     Georgia    South 152725.21 [km^2]      9468815     10006693
## 7     16       Idaho     West 216512.66 [km^2]      1526797      1616547
## 8     18     Indiana  Midwest  93648.40 [km^2]      6417398      6568645
## 9     20      Kansas  Midwest 213037.08 [km^2]      2809329      2892987
## 10    22   Louisiana    South 122345.76 [km^2]      4429940      4625253
##                          geometry
## 1  MULTIPOLYGON (((-88.20006 3...
## 2  MULTIPOLYGON (((-114.7196 3...
## 3  MULTIPOLYGON (((-109.0501 4...
## 4  MULTIPOLYGON (((-73.48731 4...
## 5  MULTIPOLYGON (((-81.81169 2...
## 6  MULTIPOLYGON (((-85.60516 3...
## 7  MULTIPOLYGON (((-116.916 45...
## 8  MULTIPOLYGON (((-87.52404 4...
## 9  MULTIPOLYGON (((-102.0517 4...
## 10 MULTIPOLYGON (((-92.01783 2...
31 / 73
data("us_states") # from spData package
us_states
## Simple feature collection with 49 features and 6 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -124.7042 ymin: 24.55868 xmax: -66.9824 ymax: 49.38436
## Geodetic CRS:  NAD83
## First 10 features:
##    GEOID        NAME   REGION             AREA total_pop_10 total_pop_15
## 1     01     Alabama    South 133709.27 [km^2]      4712651      4830620
## 2     04     Arizona     West 295281.25 [km^2]      6246816      6641928
## 3     08    Colorado     West 269573.06 [km^2]      4887061      5278906
## 4     09 Connecticut Norteast  12976.59 [km^2]      3545837      3593222
## 5     12     Florida    South 151052.01 [km^2]     18511620     19645772
## 6     13     Georgia    South 152725.21 [km^2]      9468815     10006693
## 7     16       Idaho     West 216512.66 [km^2]      1526797      1616547
## 8     18     Indiana  Midwest  93648.40 [km^2]      6417398      6568645
## 9     20      Kansas  Midwest 213037.08 [km^2]      2809329      2892987
## 10    22   Louisiana    South 122345.76 [km^2]      4429940      4625253
##                          geometry
## 1  MULTIPOLYGON (((-88.20006 3...
## 2  MULTIPOLYGON (((-114.7196 3...
## 3  MULTIPOLYGON (((-109.0501 4...
## 4  MULTIPOLYGON (((-73.48731 4...
## 5  MULTIPOLYGON (((-81.81169 2...
## 6  MULTIPOLYGON (((-85.60516 3...
## 7  MULTIPOLYGON (((-116.916 45...
## 8  MULTIPOLYGON (((-87.52404 4...
## 9  MULTIPOLYGON (((-102.0517 4...
## 10 MULTIPOLYGON (((-92.01783 2...
32 / 73

Plotting roads and other information

Click here

33 / 73

Map of Australia

library(ozmaps)
library(sf)
oz_states <- ozmaps::ozmap_states
oz_states
## Simple feature collection with 9 features and 1 field
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 105.5507 ymin: -43.63203 xmax: 167.9969 ymax: -9.229287
## Geodetic CRS:  GDA94
## # A tibble: 9 × 2
##   NAME                                                                  geometry
## * <chr>                                                       <MULTIPOLYGON [°]>
## 1 New South Wales              (((150.7016 -35.12286, 150.6611 -35.11782, 150.6…
## 2 Victoria                     (((146.6196 -38.70196, 146.6721 -38.70259, 146.6…
## 3 Queensland                   (((148.8473 -20.3457, 148.8722 -20.37575, 148.85…
## 4 South Australia              (((137.3481 -34.48242, 137.3749 -34.46885, 137.3…
## 5 Western Australia            (((126.3868 -14.01168, 126.3625 -13.98264, 126.3…
## 6 Tasmania                     (((147.8397 -40.29844, 147.8902 -40.30258, 147.8…
## 7 Northern Territory           (((136.3669 -13.84237, 136.3339 -13.83922, 136.3…
## 8 Australian Capital Territory (((149.2317 -35.222, 149.2346 -35.24047, 149.271…
## 9 Other Territories            (((167.9333 -29.05421, 167.9188 -29.0344, 167.93…
34 / 73

Obtain the following map

35 / 73
ggplot(oz_states) + 
  geom_sf() + 
  coord_sf()

36 / 73

coord_sf()

  • ensures that every layer in the plot uses the same projection
37 / 73

Mark following locations on the map

oz_capitals <- tibble::tribble( 
  ~city,           ~lat,     ~lon,
  "Sydney",    -33.8688, 151.2093,  
  "Melbourne", -37.8136, 144.9631, 
  "Brisbane",  -27.4698, 153.0251, 
  "Adelaide",  -34.9285, 138.6007, 
  "Perth",     -31.9505, 115.8605, 
  "Hobart",    -42.8821, 147.3272, 
  "Canberra",  -35.2809, 149.1300, 
  "Darwin",    -12.4634, 130.8456, 
)
oz_capitals
## # A tibble: 8 × 3
##   city        lat   lon
##   <chr>     <dbl> <dbl>
## 1 Sydney    -33.9  151.
## 2 Melbourne -37.8  145.
## 3 Brisbane  -27.5  153.
## 4 Adelaide  -34.9  139.
## 5 Perth     -32.0  116.
## 6 Hobart    -42.9  147.
## 7 Canberra  -35.3  149.
## 8 Darwin    -12.5  131.
38 / 73

Add other geoms

ggplot() + 
  geom_sf(data = oz_states, colour = "black", fill = NA) + 
  geom_point(data = oz_capitals, mapping = aes(x = lon, y = lat), colour = "red")

39 / 73

Add other geoms

ggplot() + 
  geom_sf(data = oz_states, colour = "black", fill = NA) + 
  geom_point(data = oz_capitals, mapping = aes(x = lon, y = lat), colour = "red") + coord_sf()

40 / 73

Your turn: Draw the electoral boundaries

oz_votes <- rmapshaper::ms_simplify(ozmaps::abs_ced)
41 / 73
library(ggplot2)
library(sf)
nc <- sf::st_read(system.file("shape/nc.shp", package = "sf"), quiet = TRUE)
point <- data.frame(id = "hi", lat = 36, long = -80) # point inside NC
nc
## Simple feature collection with 100 features and 14 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
## Geodetic CRS:  NAD27
## First 10 features:
##     AREA PERIMETER CNTY_ CNTY_ID        NAME  FIPS FIPSNO CRESS_ID BIR74 SID74
## 1  0.114     1.442  1825    1825        Ashe 37009  37009        5  1091     1
## 2  0.061     1.231  1827    1827   Alleghany 37005  37005        3   487     0
## 3  0.143     1.630  1828    1828       Surry 37171  37171       86  3188     5
## 4  0.070     2.968  1831    1831   Currituck 37053  37053       27   508     1
## 5  0.153     2.206  1832    1832 Northampton 37131  37131       66  1421     9
## 6  0.097     1.670  1833    1833    Hertford 37091  37091       46  1452     7
## 7  0.062     1.547  1834    1834      Camden 37029  37029       15   286     0
## 8  0.091     1.284  1835    1835       Gates 37073  37073       37   420     0
## 9  0.118     1.421  1836    1836      Warren 37185  37185       93   968     4
## 10 0.124     1.428  1837    1837      Stokes 37169  37169       85  1612     1
##    NWBIR74 BIR79 SID79 NWBIR79                       geometry
## 1       10  1364     0      19 MULTIPOLYGON (((-81.47276 3...
## 2       10   542     3      12 MULTIPOLYGON (((-81.23989 3...
## 3      208  3616     6     260 MULTIPOLYGON (((-80.45634 3...
## 4      123   830     2     145 MULTIPOLYGON (((-76.00897 3...
## 5     1066  1606     3    1197 MULTIPOLYGON (((-77.21767 3...
## 6      954  1838     5    1237 MULTIPOLYGON (((-76.74506 3...
## 7      115   350     2     139 MULTIPOLYGON (((-76.00897 3...
## 8      254   594     2     371 MULTIPOLYGON (((-76.56251 3...
## 9      748  1190     2     844 MULTIPOLYGON (((-78.30876 3...
## 10     160  2038     5     176 MULTIPOLYGON (((-80.02567 3...
42 / 73

Data

ggplot() +
  geom_sf(data = nc) +
  coord_sf(crs = st_crs(3347))

Points

ggplot() +
  geom_sf(data = nc) +
  geom_point(data = point, aes(x = long, y = lat))

43 / 73

Data and Points

ggplot() +
  geom_sf(data = nc) +
  geom_point(data = point, aes(x = long, y = lat)) +
  coord_sf(crs = st_crs(3347))

44 / 73

NAD83 / Statistics Canada Lambert

Coordinate Reference System Formats: Three common formats include.

  • proj.4

  • EPSG

  • Well-known Text (WKT) formats.

45 / 73
point <- st_as_sf(point, coords = c("long", "lat"), crs = 4326)
ggplot() +
  geom_sf(data = nc) +
  geom_sf(data = point) +
  coord_sf(crs = st_crs(3347))

coord_sf knows what to transform from

46 / 73

coord_sf can convert layers to common projections but they need to be sf objects for it to know how to do so.

47 / 73

Projecting your data

library(sf)
library(ggplot2)
# devtools::install_github("hrbrmstr/albersusa")
library(albersusa)
## Please note that 'maptools' will be retired during October 2023,
## plan transition at your earliest convenience (see
## https://r-spatial.org/r/2023/05/15/evolution4.html and earlier blogs
## for guidance);some functionality will be moved to 'sp'.
##  Checking rgeos availability: TRUE
## rgeos version: 0.6-4, (SVN revision 699)
##  GEOS runtime version: 3.10.2-CAPI-1.16.0 
##  Please note that rgeos will be retired during October 2023,
## plan transition to sf or terra functions using GEOS at your earliest convenience.
## See https://r-spatial.org/r/2023/05/15/evolution4.html for details.
##  GEOS using OverlayNG
##  Linking to sp version: 2.0-0 
##  Polygon checking: TRUE
48 / 73
crs_use = "+proj=laea +lat_0=30 +lon_0=-95"
d_points = data.frame(long = c(-110, -103, -84), 
                      lat  = c(45, 40, 41))
49 / 73
A <- ggplot(data = usa_sf()) +
  geom_sf() + 
  geom_point(data = d_points, 
             aes(x = long, y = lat), 
             color = "red", size = 5) + 
  theme_minimal() + 
  ggtitle("(A) right point position, wrong projection")
50 / 73
A
## 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()

51 / 73
B <- ggplot(data = usa_sf()) +
  geom_sf() + 
  geom_point(data = d_points, 
             aes(x = long, y = lat), 
             color = "red", size = 5) + 
  coord_sf(crs = crs_use) + 
  theme_minimal() + 
  ggtitle("(B) right projection, wrong points using geom_point()")
52 / 73
B
## 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()

53 / 73
C <- ggplot() + 
  geom_sf(data = usa_sf()) + 
  geom_sf(data = st_as_sf(d_points,
                          coords = c("long", "lat"), crs = crs_use), 
          color = "red", size = 5) + 
  coord_sf(crs = crs_use) + 
  theme_minimal() + 
  ggtitle("(C) right projection, wrong points using geom_sf() points")
54 / 73
C
## old-style crs object detected; please recreate object with a recent sf::st_crs()

55 / 73
#remotes::install_github("ropensci/USAboundaries")
#install.packages("USAboundariesData", repos = "https://ropensci.r-universe.dev", type = "source")
library(USAboundaries)
library(tidyverse)
crs_use <- "+proj=laea +lat_0=30 +lon_0=-95"
usa_sf <- us_boundaries(type="state", resolution = "low") %>% 
  dplyr::filter(!state_abbr %in% c("PR", "AK", "HI")) %>% 
  st_transform(crs = crs_use)
d_points <- data.frame(long = c(-110, -103, -84), 
                      lat  = c(45, 40, 41)) %>% 
  st_as_sf(coords = c("long", "lat"), crs = 4326) %>% 
  st_transform(crs = crs_use)
D <- ggplot(data = usa_sf) +
  geom_sf() + 
  geom_sf(data = d_points, 
             color = "red", size = 5) + 
  theme_minimal()
56 / 73
D
## function (expr, name) 
## .External(C_doD, expr, name)
## <bytecode: 0x7fbf8041fc40>
## <environment: namespace:stats>
57 / 73

st_as_sf(): what coordinate system do the coordinates you just gave me use?

st_transform(): what coordinate system should I transform those coordinates into?

Analogy: if you were translating a book,

the crs in st_as_sf() is "what language is the book written in?"

st_transform() is "what language am I translating this to?"

58 / 73

What does CRS 4326 mean?

  • The Three Most Common Coordinate Reference Systems

  • WGS84 and EPSG:4326 are geographic CRS.

  • It means that they use latitude and longitude coordinates to specify a location on the surface of the earth.

59 / 73

Projecting your data

Additional reading

Click here

60 / 73

Interactive chart

library(tidyverse)
library(sf)
library(viridis)
#devtools::install_github("tylermorganwall/rayshader")
library(rayshader)
sightings = st_read(dsn = "~/Dropbox/R Code Repository/Sasquatch/WA_OR_Bigfoot.shp")
sightings = mutate(sightings, Sightings = Join_Count)
outputFile <- file.path("~/Dropbox/R Code Repository/Sasquatch/", "sasquatch_sightings.mp4")

Resources: https://geospatialtraining.com/3d-mapping-of-sasquatch-sightings-in-r-with-rayshader/

61 / 73
sightings
62 / 73
ggSightings <- ggplot(data = sightings) + 
  geom_sf(aes(fill = Sightings)) + 
  scale_fill_viridis("Sightings") +
  ggtitle("Sasquatch Sightings - Washington and Oregon") +
  theme_bw()
63 / 73
ggSightings
64 / 73
plot_gg(ggSightings,multicore=TRUE,width=5,height=5,scale=250,windowsize=c(1400,866),
        zoom = 0.55, phi = 30)
render_snapshot()
render_movie(filename = outputFile)
65 / 73

2. Polygon map

66 / 73
# devtools::install_github("UrbanInstitute/urbnmapr")
library(tidyverse)
library(magrittr)
## 
## Attaching package: 'magrittr'
## The following object is masked from 'package:raster':
## 
##     extract
## The following object is masked from 'package:purrr':
## 
##     set_names
## The following object is masked from 'package:tidyr':
## 
##     extract
library(urbnmapr)
library(usmap) # to obtain population data
library(ggthemes)
library(scales)
## 
## Attaching package: 'scales'
## The following object is masked from 'package:viridis':
## 
##     viridis_pal
## The following object is masked from 'package:purrr':
## 
##     discard
## The following object is masked from 'package:readr':
## 
##     col_factor
library(geofacet)
67 / 73
statepop1 <- statepop %>% rename(state_name = full)
statepop2 <- full_join(statepop1, states, by = "state_name")
head(statepop2)
## # A tibble: 6 × 12
##   fips  abbr  state_name pop_2015  long   lat order hole  piece group state_fips
##   <chr> <chr> <chr>         <dbl> <dbl> <dbl> <int> <lgl> <fct> <fct> <chr>     
## 1 01    AL    Alabama     4858979 -88.5  31.9     1 FALSE 1     01.1  01        
## 2 01    AL    Alabama     4858979 -88.5  31.9     2 FALSE 1     01.1  01        
## 3 01    AL    Alabama     4858979 -88.5  31.9     3 FALSE 1     01.1  01        
## 4 01    AL    Alabama     4858979 -88.5  32.0     4 FALSE 1     01.1  01        
## 5 01    AL    Alabama     4858979 -88.5  32.0     5 FALSE 1     01.1  01        
## 6 01    AL    Alabama     4858979 -88.5  32.1     6 FALSE 1     01.1  01        
## # ℹ 1 more variable: state_abbv <chr>
68 / 73
ggplot() + 
  geom_point(data = statepop2, mapping = aes(x = long, y = lat, group = group), fill = "grey", color = "black")

69 / 73
ggplot() + 
  geom_polygon(data = statepop2, mapping = aes(x = long, y = lat, group = group),fill = "grey", color = "black")

70 / 73

Visit here: https://tstdataviz.netlify.app/slides/spatiotemporal.html

71 / 73

Moving forward

Visualizing using webmaps

library(leaflet)

Visualization of raster data

library(raster)
72 / 73
73 / 73

R packages for spatial data analysis

  • rgdal

  • sp

  • rgeos

  • raster

  • sf

  • tmap

  • leaflet

  • ggmap

  • maptools

  • gstat

  • spatstat

  • stars

  • geosphere

  • RgoogleMaps

  • rasterVis

2 / 73
Paused

Help

Keyboard shortcuts

↑, ←, Pg Up, k Go to previous slide
↓, →, Pg Dn, Space, j Go to next slide
Home Go to first slide
End Go to last slide
Number + Return Go to specific slide
b / m / f Toggle blackout / mirrored / fullscreen mode
c Clone slideshow
p Toggle presenter mode
s Start & Stop the presentation timer
t Reset the presentation timer
?, h Toggle this help
Esc Back to slideshow