Practice 1

  1. Download the coastline shapefile on your own computer.
  2. Read the shapefile with the st_read() function.

If you want, you can try with your own data or download other shapefiles from https://www.naturalearthdata.com/

Solution

# if sf is not installed
# install.packages("sf")
library(sf)

# download the data 
download.file("http://www.naturalearthdata.com/http//www.naturalearthdata.com/download/10m/physical/ne_10m_coastline.zip", 
              destfile = 'assets/data/coastlines.zip')

# unzip the file
unzip(zipfile = "assets/data/coastlines.zip", 
      exdir = 'assets/data/ne-coastlines-10m')

# WARNING: You have to define you own path in destfile argument and exdir.

# Read the file
coast <- read_sf("assets/data/ne-coastlines-10m")

Make a visual check of the spatial R object

# if mapview is not installed
# install.packages("mapview")
library(mapview)
mapview(coast)

Hint: Manage your active directory

To import files on R from your computer, you have to set properly active working directory: where R is reading and writing files on your computer. Some hints / functions to deal with that aspect.

getwd() # Print the path of the active folder / working directory
setwd(choose.dir()) 
dir() # function which list all the files in the working directory.
  1. setwd() is a function to set a new active directory
  2. choose.dir() works exclusively on windows, This function opens a dialog box to set your custom working directory

Practice 2

  1. Download and read with sf the canadian cities CSV file.
  2. Select all cities from your province
  3. Select the 10 most populated cities in your province

Solution

library(dplyr)
# 1. Read file
ca_cities <- read.csv("https://simplemaps.com/static/data/country-cities/ca/ca.csv")
ca_cities <- st_as_sf(ca_cities, coords = c("lng", "lat"), crs = 4326)

# 2. Select all cities from your province
head(ca_cities)
## Simple feature collection with 6 features and 7 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: -123.1333 ymin: 43.66667 xmax: -73.58333 ymax: 53.55
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
##        city country iso2            admin capital population
## 1   Toronto  Canada   CA          Ontario   admin    5213000
## 2  Montréal  Canada   CA           Québec            3678000
## 3 Vancouver  Canada   CA British Columbia            2313328
## 4    Ottawa  Canada   CA          Ontario primary    1145000
## 5   Calgary  Canada   CA          Alberta            1110000
## 6  Edmonton  Canada   CA          Alberta   admin    1058000
##   population_proper                   geometry
## 1           3934421 POINT (-79.41667 43.66667)
## 2           2356556     POINT (-73.58333 45.5)
## 3            603502    POINT (-123.1333 49.25)
## 4            812129     POINT (-75.7 45.41667)
## 5            915322 POINT (-114.0833 51.08333)
## 6            712391       POINT (-113.5 53.55)
qc_cities <- filter(ca_cities, admin == "Québec")
mapview(qc_cities)
# 3. Sort the cities dataframe by population
ca_cites <- ca_cities[order(ca_cities$population),]
# Print the 10 first cities
ca_cities[1:10,]
## Simple feature collection with 10 features and 7 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: -123.1333 ymin: 43.2561 xmax: -71.25 ymax: 53.55
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
##         city country iso2            admin capital population
## 1    Toronto  Canada   CA          Ontario   admin    5213000
## 2   Montréal  Canada   CA           Québec            3678000
## 3  Vancouver  Canada   CA British Columbia            2313328
## 4     Ottawa  Canada   CA          Ontario primary    1145000
## 5    Calgary  Canada   CA          Alberta            1110000
## 6   Edmonton  Canada   CA          Alberta   admin    1058000
## 7   Hamilton  Canada   CA          Ontario             721053
## 8   Winnipeg  Canada   CA         Manitoba   admin     632063
## 9     Québec  Canada   CA           Québec   admin     624177
## 10    Oshawa  Canada   CA          Ontario             450963
##    population_proper                   geometry
## 1            3934421 POINT (-79.41667 43.66667)
## 2            2356556     POINT (-73.58333 45.5)
## 3             603502    POINT (-123.1333 49.25)
## 4             812129     POINT (-75.7 45.41667)
## 5             915322 POINT (-114.0833 51.08333)
## 6             712391       POINT (-113.5 53.55)
## 7             519949  POINT (-79.85748 43.2561)
## 8             575313 POINT (-97.16667 49.88333)
## 9             528595        POINT (-71.25 46.8)
## 10            247989     POINT (-78.86667 43.9)

Practice 3

  1. Read with sf the canadian cities CSV file.
  2. Import with getData the canadian digital elevation model (elevation raster)
  3. Project the elevation raster and the canadian cities in NAD83 (raster::projectRaster() & sf::st_transform())
  4. Make a buffer of 10 kms around each canadian cities (sf::st_buffer())
  5. Compute the elevation mean per buffer (raster::extract())

Solution

library(raster)
## Loading required package: sp
## 
## Attaching package: 'raster'
## The following object is masked from 'package:dplyr':
## 
##     select
# 1. Read file
ca_cities <- read.csv("https://simplemaps.com/static/data/country-cities/ca/ca.csv")
ca_cities <- st_as_sf(ca_cities, coords = c("lng", "lat"), crs = 4326)

# 2. Import the digital elevation model
altCAN <- getData(name = "alt", country = "CAN", path = "assets/data") # Coarse resolution

# 3 Print the current CRS of elevation raster 
proj4string(altCAN)
## [1] "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"
# Print the current CRS of the canadian cities
st_crs(ca_cities)
## Coordinate Reference System:
##   EPSG: 4326 
##   proj4string: "+proj=longlat +datum=WGS84 +no_defs"
# Reproject ca_cities and altCan in 3347
# https://spatialreference.org/ref/epsg/3347/
# With sf
ca_cities <- st_transform(ca_cities, 3347)
# With raster
# This takes a while
altCAN <- projectRaster(altCAN, crs="+proj=lcc +lat_1=49 +lat_2=77 +lat_0=63.390675 +lon_0=-91.86666666666666 +x_0=6200000 +y_0=3000000 +ellps=GRS80 +datum=NAD83 +units=m +no_defs;")

# 4. Generate the buffers
cities_buf <- st_buffer(ca_cities, 10000)
mapview(cities_buf)
# 5. Crop elevation raster with buffers and compute the mean elevation
mu_elev_cities <- extract(altCAN, cities_buf, fun = mean, na.rm = TRUE)
head(mu_elev_cities)
##            [,1]
## [1,]  115.11302
## [2,]   32.92916
## [3,]   39.87088
## [4,]   73.54229
## [5,] 1105.83312
## [6,]  667.81916
length(mu_elev_cities)
## [1] 247
# Visual check of the result
ca_cities$mu_elev_cities <- as.numeric(mu_elev_cities)
mapview(ca_cities, zcol = "mu_elev_cities", cex = "mu_elev_cities")