class: title-slide, middle <img src="assets/img/pr.png" width="80px" style="padding-right:10px; border-right: 2px solid #FAFAFA;"></img> <img src="assets/img/UdeS_blanc.png" width="250px" ></img> # Handling spatial data with R .instructors[ GSFE01 - F. Guillaume Blanchet & Steve Vissault ] --- class: clear, middle ##
Slides available at: https://steveviss.github.io/PR-GSFE01/basics ##
Solutions of all practices: https://steveviss.github.io/PR-GSFE01/basics/practice.html --- class: inverse, center, middle # Review basic R instructions <html><div style='float:left'></div><hr color='#EB811B' size=1px width=720px></html> --- # R basics ## Object assignation Code in the black area ```r x <- 3 x ``` Output starting with `##` ``` ## [1] 3 ``` `x` is an object whichin we assign (`<-`) a numerical value of `3` --- # R basics R can be used as a calculator ```r 1.2 + -0.654 * 32 ``` ``` ## [1] -19.728 ``` Values stored in the R environment (through objects) can be reused ```r y <- x + -0.654 * 32 y ``` ``` ## [1] -17.928 ``` --- # Logical operators ### Several options ``` x == y # equal x != y # not equal x < y # smaller than x >= y # greater than or equal ``` ### Example ```r x == y ``` ``` ## [1] FALSE ``` ```r # Return logical values TRUE / FALSE ``` --- class: inverse, center, middle # Vector <html><div style='float:left'></div><hr color='#EB811B' size=1px width=720px></html> --- # Vector ### Vector declaration ```r x <- c("a", "b", "c") x ``` ``` ## [1] "a" "b" "c" ``` ```r # Other way with numerical values 1:3 ``` ``` ## [1] 1 2 3 ``` -- ### Can we mixed type of values (ex. numeric with text)? ```r x <- c("a",1,2) ``` -- Yes, but the entire vector will be convert as a character string vector --- # Vector position ### Accessing specific position in vector ```r x <- c("a", "b", "c") x[1] ``` ``` ## [1] "a" ``` ```r x[c(2,1)] ``` ``` ## [1] "b" "a" ``` --- # Vector manipulations ### Assigning or removing items in the array ```r x[1] <- "z" x ``` ``` ## [1] "z" "b" "c" ``` ```r x[-c(2,1)] ``` ``` ## [1] "c" ``` ### Vector algebra ```r 3 * c(11,22,33) ``` ``` ## [1] 33 66 99 ``` --- # Vector Special values ```r c(NA, 1/0, 0/0) ``` ``` ## [1] NA Inf NaN ``` --- class: inverse, center, middle # Matrix <html><div style='float:left'></div><hr color='#EB811B' size=1px width=720px></html> --- # Matrix ### Matrix declaration ```r m <- matrix(20:23,nrow=2,ncol=2) m ``` ``` ## [,1] [,2] ## [1,] 20 22 ## [2,] 21 23 ``` ### Accessing specific positions in matrix ```r m[1,] m[,2] m[1,2] ``` ``` ## [1] 20 22 ## [1] 22 23 ## [1] 22 ``` --- # Matrix algebra ```r m <- matrix(c(1,1,2,2,3,3),nrow=3,ncol=2) m ``` ``` ## [,1] [,2] ## [1,] 1 2 ## [2,] 1 3 ## [3,] 2 3 ``` ```r m %*% t(m) ``` ``` ## [,1] [,2] [,3] ## [1,] 5 7 8 ## [2,] 7 10 11 ## [3,] 8 11 13 ``` --- # Different columns type ### Can we have two columns type within a matrix? ```r matrix(c("a","b",2, 1),nrow=2,ncol=2) ``` ``` ## [,1] [,2] ## [1,] "a" "2" ## [2,] "b" "1" ``` Same as for a vector, the entire matrix will be coerce into character string. -- .pull-left[ So, how can we get a table with two columns with different types of information? ] -- .pull-right[ <img src="./assets/img/data_structures2.png" width="60%" style="display: block; margin: auto;" /> ] --- class: inverse, center, middle # Dataframe <html><div style='float:left'></div><hr color='#EB811B' size=1px width=720px></html> --- # Dataframe Table generated by combining several vectors of any type (text, date, numeric etc.) ```r df <- data.frame( col_text = c("cat","dog","parrot","guinea pig"), col_date = as.Date(c("2019-11-12","2019-11-11","2019-11-10","2019-11-09")), col_numerical = c(11,22,33,44), stringsAsFactors = FALSE ) ``` ```r df ``` ``` ## col_text col_date col_numerical ## 1 cat 2019-11-12 11 ## 2 dog 2019-11-11 22 ## 3 parrot 2019-11-10 33 ## 4 guinea pig 2019-11-09 44 ``` --- # Dataframe ### Accessing specific element ```r # Row df[1,] ``` ``` ## col_text col_date col_numerical ## 1 cat 2019-11-12 11 ``` ```r # Column df[,1] ``` ``` ## [1] "cat" "dog" "parrot" "guinea pig" ``` ```r # Position df[1,1] ``` ``` ## [1] "cat" ``` --- # Dataframe ### Accessing specific columns ```r df[,c("col_text","col_date")] ``` ``` ## col_text col_date ## 1 cat 2019-11-12 ## 2 dog 2019-11-11 ## 3 parrot 2019-11-10 ## 4 guinea pig 2019-11-09 ``` ```r df$col_date ``` ``` ## [1] "2019-11-12" "2019-11-11" "2019-11-10" "2019-11-09" ``` -- With dataframe, all columns have to have the same length. Want more flexibility : `list()`. --- class: inverse, center, middle # List <html><div style='float:left'></div><hr color='#EB811B' size=1px width=720px></html> --- # List A list store collections of objects (single value, vector or matrix). ```r l <- list( v = c(1,2,3,4), m = matrix(c(1,1,1,1),nrow=2,ncol=2) ) ``` ```r l ``` ``` ## $v ## [1] 1 2 3 4 ## ## $m ## [,1] [,2] ## [1,] 1 1 ## [2,] 1 1 ``` --- # List ```r l <- list( v = c(1,2,3,4), m = matrix(c(1,1,1,1),nrow=2,ncol=2) ) ``` ```r l$v #equivalent: l[['v']] ``` ``` ## [1] 1 2 3 4 ``` --- # List ```r l <- list( v = c(1,2,3,4), m = matrix(c(1,1,1,1),nrow=2,ncol=2) ) ``` ```r l$m #equivalent: l[['m']] ``` ``` ## [,1] [,2] ## [1,] 1 1 ## [2,] 1 1 ``` --- # List ```r l <- list( v = c(1,2,3,4), m = matrix(c(1,1,1,1),nrow=2,ncol=2) ) ``` ```r l$v[2] ``` ``` ## [1] 2 ``` --- class: inverse, center, middle # Wrap up <html><div style='float:left'></div><hr color='#EB811B' size=1px width=720px></html> --- # Wrap up ### Different object types <img src="./assets/img/data_structures.png" width="60%" style="display: block; margin: auto;" /> .font70[ From https://mgimond.github.io/ES218/Week02a.html ] --- # R basics ### Functions .pull-left[ ```r Add <- function(x, y) { return(x + y) } Add(10, 3) ``` ``` ## [1] 13 ``` ] .pull-right[ <img src="https://qcbsrworkshops.github.io/workshop01/workshop01-en/images/function_V2_en.png" width="40%" style="display: block; margin: auto;" /> ] --- class: clear, middle # R Universe R contains a **plethora of packages**. Each package expose a set of functions. <img src="./assets/img/CRAN.png" width="80%" style="display: block; margin: auto;" /> From https://www.rdocumentation.org/trends --- class: clear, middle # Interact with the R universe ```r # Install packages install.packages("sf") # Then access to a specific function sf::st_crop() # Or load the entire package library(sf) ``` --- class: clear, middle # Navigate in the R universe Search in the documentation https://www.rdocumentation.org/ or https://rdrr.io/ Search or get help on specific functions ```r ?plot # Search within your loaded packages ??raster # Search among packages (loaded and unloaded) ``` --- class: clear, middle, center # R Universe <img src="./assets/img/cloud-packages.png" width="65%" style="display: block; margin: auto;" /> --- class: clear, middle, center # Spatial galaxy <img src="./assets/img/cloud-packages-zoomin.png" width="65%" style="display: block; margin: auto;" /> --- class: inverse, center, middle # Spatial galaxy <html><div style='float:left'></div><hr color='#EB811B' size=1px width=720px></html> --- class: clear, middle # What's coming next? 1. Overview of spatial objects 2. Manipulate spatial vectors (POINT, LINE, POLYGON) 4. Spatial operations on vectors 5. Manipulate rasters 6. Spatial operations on rasters --- class: inverse, center, middle # Overview of spatial objects <html><div style='float:left'></div><hr color='#EB811B' size=1px width=720px></html> --- class: middle, clear # Spatial vector <img src="./assets/img/vectors.png" width="100%" style="display: block; margin: auto;" /> --- class: middle, clear # Rasters <img src="./assets/img/raster.png" width="100%" style="display: block; margin: auto;" /> --- class: inverse, center, middle # Coordinate Reference System (CRS) <html><div style='float:left'></div><hr color='#EB811B' size=1px width=720px></html> --- # Coordinate Reference System (CRS) Projections are defined as coordinate reference system (CRS) ### Geographic (or unprojected) CRS <img src="./assets/img/projection1.png" width="60%" style="display: block; margin: auto;" /> ### Projected CRS <img src="./assets/img/projection2.png" width="60%" style="display: block; margin: auto;" /> --- # Coordinate Reference System (CRS) ### Geographic (or unprojected) CRS - Latitude and longitude, i.e. angles measured from the Earth’s center to a point on the Earth’s surface - 3-D representation of Earth (sphere or ellipsoid) - Distance in geographic CRSs are therefore measured in degrees, not meters - Lon/Lat locate any points on Earth’s surface, but are not uniform units of measure ### Projected CRS - Uses Cartesian coordinates, Easting and Northing (x and y) **typically in meters** - 2-D representation of Earth - All projected CRSs are based on a geographic CRS - Different mathematical formulas (i.e. projections) can transform the 3-D globe to a 2-D map --- # Geographic vs projected CRS <img src="./assets/img/proj_vs_unproj.png" width="80%" style="display: block; margin: auto;" /> The representation may not be very different at fine spatial scale, **but it is very different at broader spatial scales.** --- # Coordinate Reference System (CRS) - Each CRS has a spatial reference ID called SRID or EPSG - In R, the notation used to describe the CRS is **proj4string** from the PROJ.4 library. It looks like this: `+init=epsg:4326 +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs` Knowing the EPSG/SRID of the projection can thus become handy. Luckily, there are very good resources on the web for this. To search for a specific CRS, have a look at [spatialreference.org](https://spatialreference.org/) exemple: https://spatialreference.org/ref/epsg/4326/ [Overview of Coordinate Reference Systems (CRS)](https://www.nceas.ucsb.edu/~frazier/RSpatialGuides/OverviewCoordinateReferenceSystems.pdf) --- # Coordinate Reference System (CRS) There are also R functions designed to use this information to convert EPSG/SRID for R objects - `rgdal::CRS("+init=epsg:4326")` - `sp::CRS("+init=epsg:4326")` - `sf::st_crs(4326)` --- class: inverse, center, middle # Manipulate spatial vectors <html><div style='float:left'></div><hr color='#EB811B' size=1px width=720px></html> Points, Lines & Polygons --- class: middle, clear # R packages for spatial vectors - `sp` package provides classes and methods for dealing with spatial data. - `maptools` & `rgeos` packages give access to functions to manipulate and analyze spatial vectors - `sf` is newer. It provides Simple Features for R, in compliance with the [Open Geospatial Consortium](https://www.opengeospatial.org/) (OGC) Simple Feature standard. -- In this course, we will focus mainly on `sf`. --- class: inverse, center, middle # Why `sf`? <html><div style='float:left'></div><hr color='#EB811B' size=1px width=720px></html> --- # Why sf? Why use the `sf` package when `sp` is already tried and tested? 1. Simple features refers to a formal standard (ISO 19125-1:2004) that describes how objects in the real world can be represented in computers 2. Successor of `sp` 3. `sf` incorporates the functionality of the 3 mai spatial data manipulation packages in a single, cohesive whole: - `sp` for the class system; - `rgdal` for reading and writing data; - `rgeos` for spatial operations undertaken by GEOS. So, instead of learning 3 packages, we focus only on one. --- # Why sf? - `sf` objects are easy to manipulate: spatial objects are stored as data frames, with the feature geometries stored in list-columns <img src="https://user-images.githubusercontent.com/520851/50280460-e35c1880-044c-11e9-9ed7-cc46754e49db.jpg" width="70%" style="display: block; margin: auto;" /> .font70[ Illustration from [Allison Horst](https://twitter.com/allison_horst/status/1071456081308614656) ] --- # Structure of `sf` object <img src="./assets/img/str_sf.png" width="100%" style="display: block; margin: auto;" /> --- # `sfg` - Simple feature geometry <img src="./assets/img/sf1.png" width="70%" style="display: block; margin: auto;" /> --- # `sfg` - Simple feature geometry <img src="./assets/img/sf2.png" width="70%" style="display: block; margin: auto;" /> --- class: inverse, center, middle # Creating simple feature <html><div style='float:left'></div><hr color='#EB811B' size=1px width=720px></html> --- # Points Declare spatial vertices / points (same as .shp file) ```r library(sf) # Do not forget to load the library ottawa <- st_point(c(-75.69812, 45.41117)) sherbrooke <- st_point(c(-71.89908, 45.40008)) winnipeg <- st_point(c(-97.14704, 49.8844)) calgary <- st_point(c(-114.08529, 51.05011)) vancouver <- st_point(c(-123.11934, 49.24966)) ``` Set all points in spatial column and the 4326 CRS (WGS84) ```r cities <- st_sfc( list(ottawa, sherbrooke, winnipeg, calgary, vancouver), crs = 4326 ) class(cities) ``` ``` ## [1] "sfc_POINT" "sfc" ``` --- # Declare attributes ```r # Declare attribute table (same as .dbf file) attr_table <- data.frame( N = c(1236324, 212105, 804200, 1406700, 2463431), name = c("Ottawa","Sherbrooke","Winnipeg","Calgary","Vancouver") ) attr_table ``` ``` ## N name ## 1 1236324 Ottawa ## 2 212105 Sherbrooke ## 3 804200 Winnipeg ## 4 1406700 Calgary ## 5 2463431 Vancouver ``` --- # Declare CRS ```r # Find approriate projection proj <- st_crs(4326) proj ``` ``` ## Coordinate Reference System: ## EPSG: 4326 ## proj4string: "+proj=longlat +datum=WGS84 +no_defs" ``` --- # Attach attributes table and CRS ```r # Attach spatial feature, attribute table great_cities <- st_sf(attr_table, geom = cities, crs = proj) great_cities ``` ``` ## Simple feature collection with 5 features and 2 fields ## geometry type: POINT ## dimension: XY ## bbox: xmin: -123.1193 ymin: 45.40008 xmax: -71.89908 ymax: 51.05011 ## epsg (SRID): 4326 ## proj4string: +proj=longlat +datum=WGS84 +no_defs ## N name geom ## 1 1236324 Ottawa POINT (-75.69812 45.41117) ## 2 212105 Sherbrooke POINT (-71.89908 45.40008) ## 3 804200 Winnipeg POINT (-97.14704 49.8844) ## 4 1406700 Calgary POINT (-114.0853 51.05011) ## 5 2463431 Vancouver POINT (-123.1193 49.24966) ``` --- # Visual check ```r plot(great_cities[,"N"]) ``` <img src="index_files/figure-html/unnamed-chunk-50-1.png" width="432" style="display: block; margin: auto;" /> --- # Visual check ```r library(mapview) mapview(great_cities, zoom = 1) ```
-- What about spatial lines and polygons? --- # Lines Declaring lines and polygons may become complex if you want to account for fine details. ```r # Example of line composed of 3 coordinates line <- st_linestring(rbind(c(0,0),c(1,1),c(2,1))) class(line) ``` ``` ## [1] "XY" "LINESTRING" "sfg" ``` ```r # Visualise plot(line, col = "red", lwd=2) plot(st_cast(line,"MULTIPOINT"), pch = 19, add = TRUE) ``` <img src="index_files/figure-html/unnamed-chunk-52-1.png" width="288" style="display: block; margin: auto;" /> --- # Polygons Declaring lines and polygons may become complex, when accounting for fine details. ```r outer = matrix(c(0,0,10,0,10,10,0,10,0,0),ncol=2, byrow=TRUE) hole1 = matrix(c(1,1,1,2,2,2,2,1,1,1),ncol=2, byrow=TRUE) hole2 = matrix(c(5,5,5,6,6,6,6,5,5,5),ncol=2, byrow=TRUE) pts = list(outer, hole1, hole2) poly <- st_polygon(pts) plot(poly, col = "red") ``` <img src="index_files/figure-html/unnamed-chunk-53-1.png" width="216" style="display: block; margin: auto;" /> -- This is why we usually import them from shapefiles (lines & polygons) and CSV (points). --- class: inverse, center, middle # Importing and exporting simple features <html><div style='float:left'></div><hr color='#EB811B' size=1px width=720px></html> --- # Importing simple points features from CSV ```r ca_cities <- read.csv("https://simplemaps.com/static/data/country-cities/ca/ca.csv") head(ca_cities, 4) ``` ``` ## city lat lng country iso2 admin capital ## 1 Toronto 43.66667 -79.41667 Canada CA Ontario admin ## 2 Montréal 45.50000 -73.58333 Canada CA Québec ## 3 Vancouver 49.25000 -123.13333 Canada CA British Columbia ## 4 Ottawa 45.41667 -75.70000 Canada CA Ontario primary ## population population_proper ## 1 5213000 3934421 ## 2 3678000 2356556 ## 3 2313328 603502 ## 4 1145000 812129 ``` -- How to make it spatial? We need to know 3 basics informations: - Latitude (`lat`) - Longitude (`lng`), - Projection (i.e. the CRS) Unfortunately, projection (the CRS) is often hidden somewhere in the CSV file metadata. --- # Importing simple features from CSV ```r sf_ca_cities_wgs84 <- st_as_sf(ca_cities, coords = c("lng", "lat"), crs = 4326) sf_ca_cities_wgs84 ``` ``` ## Simple feature collection with 247 features and 7 fields ## geometry type: POINT ## dimension: XY ## bbox: xmin: -139.45 ymin: 42.30165 xmax: -52.66667 ymax: 82.48333 ## epsg (SRID): 4326 ## proj4string: +proj=longlat +datum=WGS84 +no_defs ## First 10 features: ## 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) ``` And it's done. --- # Importing features from files `sf` can read a large number of spatial vector formats. This is possible because of the hardwork of the OGC organisation developing international standards. The first 50 formats ```r st_drivers(what = "vector")[1:50,1] ``` ``` ## [1] PCIDSK netCDF JP2OpenJPEG PDF ## [5] ESRI Shapefile MapInfo File UK .NTF OGR_SDTS ## [9] S57 DGN OGR_VRT REC ## [13] Memory BNA CSV NAS ## [17] GML GPX LIBKML KML ## [21] GeoJSON Interlis 1 Interlis 2 OGR_GMT ## [25] GPKG SQLite OGR_DODS ODBC ## [29] WAsP PGeo MSSQLSpatial OGR_OGDI ## [33] PostgreSQL MySQL OpenFileGDB XPlane ## [37] DXF CAD Geoconcept GeoRSS ## [41] GPSTrackMaker VFK PGDUMP OSM ## [45] GPSBabel SUA OpenAir OGR_PDS ## [49] WFS SOSI ## 210 Levels: AAIGrid ACE2 ADRG AeronavFAA AIG AirSAR AmigoCloud ... ZMap ``` --- # Importing features from ESRI shapefiles Let's download the world coastlines shapefile (available [here](https://www.naturalearthdata.com/http//www.naturalearthdata.com/download/10m/physical/ne_10m_coastline.zip)) ```r # 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. ``` --- # Importing features from ESRI shapefiles ```r dir('assets/data/ne-coastlines-10m') ``` ``` ## [1] "ne_10m_coastline.cpg" "ne_10m_coastline.dbf" ## [3] "ne_10m_coastline.prj" "ne_10m_coastline.README.html" ## [5] "ne_10m_coastline.shp" "ne_10m_coastline.shx" ## [7] "ne_10m_coastline.VERSION.txt" ``` ### Valid ESRI shapefiles - `.shp`: The main file that stores the feature geometry (required). - `.shx`: The index file that stores the index of the feature geometry (required). - `.dbf`: The dBASE table that stores the attribute information of features (required). - `.prj`: The file that stores the coordinate system information (used by ArcGIS). --- # Importing features from ESRI shapefiles ```r coast <- st_read('assets/data/ne-coastlines-10m') ``` ``` ## Reading layer `ne_10m_coastline' from data source `/home/steve/Documents/Git/PRStats/PR-GSFE01/basics/assets/data/ne-coastlines-10m' using driver `ESRI Shapefile' ## Simple feature collection with 4133 features and 3 fields ## geometry type: LINESTRING ## dimension: XY ## bbox: xmin: -180 ymin: -85.22194 xmax: 180 ymax: 83.6341 ## epsg (SRID): 4326 ## proj4string: +proj=longlat +datum=WGS84 +no_defs ``` ```r par(mar=rep(0,4)) plot(st_geometry(coast)) ``` <img src="index_files/figure-html/unnamed-chunk-59-1.png" style="display: block; margin: auto;" /> --- # Importing features from GeoDataBase Simple exemple from the Community Aquatic Monitoring Program (CAMP) for the Southern Gulf of St. Lawrence (Download from [open.canada.ca](https://open.canada.ca/data/en/dataset/c4474517-3d9b-e581-a6e2-e95273f2058e)) ```r st_layers("assets/data/camp_station_summary_eng.gdb") ``` ``` ## Driver: OpenFileGDB ## Available layers: ## layer_name geometry_type features fields ## 1 camp_station_summary_eng Point 228 10 ``` ```r stations <- st_read("assets/data/camp_station_summary_eng.gdb", layer = "camp_station_summary_eng") ``` ``` ## Reading layer `camp_station_summary_eng' from data source `/home/steve/Documents/Git/PRStats/PR-GSFE01/basics/assets/data/camp_station_summary_eng.gdb' using driver `OpenFileGDB' ## Simple feature collection with 228 features and 10 fields ## geometry type: POINT ## dimension: XY ## bbox: xmin: -65.66146 ymin: 45.62596 xmax: -61.01743 ymax: 47.81772 ## epsg (SRID): 4326 ## proj4string: +proj=longlat +datum=WGS84 +no_defs ``` --- # Importing features from GeoDataBase Simple exemple from the Community Aquatic Monitoring Program (CAMP) for the Southern Gulf of St. Lawrence (Download from [open.canada.ca](https://open.canada.ca/data/en/dataset/c4474517-3d9b-e581-a6e2-e95273f2058e))
--- # Export spatial vectors We go back to the coastline object and want to extract the St. Lawrence Gulf. First, we draw the polygon of the area of interest (using http://arthur-e.github.io/Wicket/sandbox-gmaps3.html). ```r area <- st_as_sfc("POLYGON((-71.19219092922373 51.818174659518405,-55.06426124172373 51.818174659518405,-55.06426124172373 45.47097576656452,-71.19219092922373 45.47097576656452,-71.19219092922373 51.818174659518405))") # This is a WKT representation of the polygon, often used on the web, # a quick way to define a polygon area <- st_set_crs(area, 4326) cropped_coast <- st_crop(coast,area) ``` ``` ## although coordinates are longitude/latitude, st_intersection assumes that they are planar ``` ``` ## Warning: attribute variables are assumed to be spatially constant ## throughout all geometries ``` --- # Export spatial vectors ```r # Visual representation par(mar=rep(0,4)) plot(st_geometry(cropped_coast)) ``` <img src="index_files/figure-html/unnamed-chunk-63-1.png" style="display: block; margin: auto;" /> --- # Export spatial vectors Then, we save the new cropped spatial feature to an ESRI shapefile on the computer. ```r st_write(cropped_coast, dsn = "assets/data/cropped_coast", driver = "ESRI Shapefile") ``` --- class: clear, middle #
Practice 1. Download the [coastline shapefile](https://www.naturalearthdata.com/http//www.naturalearthdata.com/download/10m/physical/ne_10m_coastline.zip) 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/ --- class: inverse, center, middle # Manipulating spatial features with `sf` <html><div style='float:left'></div><hr color='#EB811B' size=1px width=720px></html> --- # Filtering based on data attributes Select cities in Manitoba ```r mb_cities <- subset(sf_ca_cities_wgs84, admin == 'Manitoba') mapview(mb_cities) ```
--- # Select a spatial feature ```r # Select the first feature mb_cities[1,] ``` ``` ## Simple feature collection with 1 feature and 7 fields ## geometry type: POINT ## dimension: XY ## bbox: xmin: -97.16667 ymin: 49.88333 xmax: -97.16667 ymax: 49.88333 ## epsg (SRID): 4326 ## proj4string: +proj=longlat +datum=WGS84 +no_defs ## city country iso2 admin capital population population_proper ## 8 Winnipeg Canada CA Manitoba admin 632063 575313 ## geometry ## 8 POINT (-97.16667 49.88333) ``` --- # Remove spatial features ```r # Remove the first 5 features mb_cities[-c(1:5),] ``` ``` ## Simple feature collection with 15 features and 7 fields ## geometry type: POINT ## dimension: XY ## bbox: xmin: -101.8833 ymin: 50.63333 xmax: -92.08333 ymax: 58.76667 ## epsg (SRID): 4326 ## proj4string: +proj=longlat +datum=WGS84 +no_defs ## First 10 features: ## city country iso2 admin capital population ## 102 Dauphin Canada CA Manitoba 9077 ## 112 Flin Flon Canada CA Manitoba 6393 ## 117 The Pas Canada CA Manitoba 6055 ## 118 Norway House Canada CA Manitoba 6000 ## 150 Gimli Canada CA Manitoba 2623 ## 152 Nelson House Canada CA Manitoba 2500 ## 180 Gillam Canada CA Manitoba 1281 ## 189 Churchill Canada CA Manitoba 1000 ## 192 Berens River Canada CA Manitoba 892 ## 193 Shamattawa Canada CA Manitoba 870 ## population_proper geometry ## 102 8418 POINT (-100.05 51.15) ## 112 6002 POINT (-101.8833 54.76667) ## 117 3802 POINT (-101.2333 53.81667) ## 118 5000 POINT (-97.83333 53.96667) ## 150 2009 POINT (-97 50.63333) ## 152 2500 POINT (-98.85 55.8) ## 180 1281 POINT (-94.7 56.35) ## 189 923 POINT (-94.16667 58.76667) ## 192 153 POINT (-97.03333 52.36667) ## 193 870 POINT (-92.08333 55.85) ``` --- # Transform CRS Change the original projection from WGS84 (SRID 4326) to NAD83 / Statistics Canada Lambert (SRID 3347). ```r sf_ca_cities_nad83 <- st_transform(sf_ca_cities_wgs84, 3347) sf_ca_cities_nad83 ``` ``` ## Simple feature collection with 247 features and 7 fields ## geometry type: POINT ## dimension: XY ## bbox: xmin: 3740150 ymin: 725605.6 xmax: 8992280 ymax: 5213670 ## epsg (SRID): 3347 ## proj4string: +proj=lcc +lat_1=49 +lat_2=77 +lat_0=63.390675 +lon_0=-91.86666666666666 +x_0=6200000 +y_0=3000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs ## First 10 features: ## 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 (7221807 929663.7) ## 2 2356556 POINT (7630595 1244080) ## 3 603502 POINT (4016310 2004568) ## 4 812129 POINT (7471164 1190060) ## 5 915322 POINT (4686217 1926913) ## 6 712391 POINT (4814910 2168989) ## 7 519949 POINT (7194844 876773.3) ## 8 575313 POINT (5820801 1542873) ## 9 528595 POINT (7760535 1438429) ## 10 247989 POINT (7260945 964674.4) ``` --- # Create buffers around cities Create buffer of 10 kms around cities. This will transform sf_ca_cities_nad83 from POINT to POLYGONS. ```r buf_10K_cities <- st_buffer(sf_ca_cities_nad83, 10000) mapview(buf_10K_cities) ```
And so many spatial operations that we cannot cover all of them... --- # Spatial operations with `sf` <img src="./assets/img/geom_ops.png" width="85%" style="display: block; margin: auto;" /> [Download `sf` cheat sheet](https://github.com/rstudio/cheatsheets/blob/master/sf.pdf) --- class: clear, middle #
Practice 1. Download and read with `sf` the [canadian cities CSV file](https://simplemaps.com/static/data/country-cities/ca/ca.csv). 2. Select all cities from your province 3. Select the 10 most populated cities in your province --- class: inverse, center, middle # Rasters manipulation <html><div style='float:left'></div><hr color='#EB811B' size=1px width=720px></html> --- # Gridded spatial data Gridded spatial data is a special case that need to be handle with the library `raster`. <img src="./assets/img/raster_pic.png" width="90%" style="display: block; margin: auto;" /> 136 raster formats supported by GDAL library (and hence by R). --- # Diversity of formats .pull-left[ **Some examples** *Classic gridded format* ``` GTiff: GeoTIFF XYZ: ASCII Gridded XYZ ``` *Images* ``` PNG: Portable Network Graphics JPEG: JPEG JFIF ``` *Multibands (often related to satellite imaging)* ``` netCDF: Network Common Data Format HDF4: Hierarchical Data Format Release 4 ``` ] .pull-right[ <img src="./assets/img/rgb.jpeg" width="100%" style="display: block; margin: auto;" /> ] --- # Load raster from file We simply use the function `raster()` to import the file in R. ```r library(raster) ocean_bottom <- raster("assets/data/OB_LR/OB_LR.tif") ocean_bottom ``` ``` ## class : RasterLayer ## band : 1 (of 3 bands) ## dimensions : 8100, 16200, 131220000 (nrow, ncol, ncell) ## resolution : 0.02222222, 0.02222222 (x, y) ## extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax) ## crs : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 ## source : /home/steve/Documents/Git/PRStats/PR-GSFE01/basics/assets/data/OB_LR/OB_LR.tif ## names : OB_LR ## values : 0, 255 (min, max) ``` --- # Load raster from file ```r image(ocean_bottom) ``` <img src="index_files/figure-html/unnamed-chunk-74-1.png" style="display: block; margin: auto;" /> --- # Retrieving free GIS data with `getData` The function `getData` in the raster package allows to freely and quickly retrieve free GIS data such as: - **GADM** - global administrative boundaries at different level of administrative subdivision - **worldclim** - global interpolated climate data - **alt** and STRM - coarse and fine resolution elevation data - **ISO3** - 3 letter ISO codes for country names. ```r # alt90 <- getData('SRTM', lon = -73.7, lat = 45.5) # Fine resolution altCAN <- getData(name = "alt", country = "CAN", path = "assets/data") # Coarse resolution altCAN ``` ``` ## class : RasterLayer ## dimensions : 4992, 10620, 53015040 (nrow, ncol, ncell) ## resolution : 0.008333333, 0.008333333 (x, y) ## extent : -141.1, -52.6, 41.6, 83.2 (xmin, xmax, ymin, ymax) ## crs : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 ## source : /home/steve/Documents/Git/PRStats/PR-GSFE01/basics/assets/data/CAN_msk_alt.grd ## names : CAN_msk_alt ## values : -108, 5800 (min, max) ``` --- # Retrieving free GIS data with `getData` ```r plot(altCAN) ``` <img src="index_files/figure-html/unnamed-chunk-76-1.png" style="display: block; margin: auto;" /> --- # Extract spatial values at specific locations To interact with raster object, we have to convert the `sf_ca_cities_wgs84` containing the canadian cities. ```r sf_ca_cities_wgs84$elev <- extract(altCAN,sf_ca_cities_wgs84) # return a vector of values with we can be attach to the sf object sf_ca_cities_wgs84$elev ``` ``` ## [1] 121 56 74 70 1081 663 92 233 74 106 332 55 242 186 ## [15] 53 484 229 578 292 56 185 185 NA 354 16 92 184 28 ## [29] NA NA 197 83 855 910 386 589 697 88 90 112 13 197 ## [43] 88 58 58 93 9 133 650 393 22 219 65 308 447 35 ## [57] 31 440 547 392 94 170 NA 297 632 311 247 69 306 44 ## [71] 63 526 243 693 932 153 88 741 123 503 741 NA 629 759 ## [85] 476 208 123 100 759 840 212 338 706 299 36 1055 226 264 ## [99] 569 17 214 293 382 471 5 1397 24 271 204 7 10 337 ## [113] 397 502 NA 23 267 222 478 635 807 360 175 245 563 300 ## [127] 271 680 6 550 18 23 17 69 1060 157 10 361 315 3 ## [141] 405 132 245 38 19 565 202 430 696 223 550 242 1 46 ## [155] 642 208 288 369 253 NA 3 372 10 NA NA 12 191 8 ## [169] NA 39 NA 172 NA 324 5 1 337 NA 336 137 NA 1551 ## [183] 193 NA NA 15 NA 17 7 1 NA 222 76 4 1 67 ## [197] 693 NA NA 285 NA 98 24 4 40 21 156 209 NA NA ## [211] 337 518 47 61 202 9 272 2 60 NA 3 993 11 145 ## [225] 338 396 122 7 8 37 87 187 52 NA 239 38 1 NA ## [239] 241 109 156 299 787 NA 253 227 311 ``` --- # Crop and mask rasters `crop()` and `mask()` reduce object memory use and therefore computational time for raster manipulations. `crop()` will decrease the extent of a raster, `mask()` will set to NA values outside a polygon boundary. ```r # Get Quebec shapefile can <- raster::getData("GADM", country = "CAN", level = 1, path = "assets/data") # convert to sf can <- st_as_sf(can) # Subset Quebec qc <- subset(can, NAME_1 == "Québec") alt_qc <- crop(altCAN, qc) alt_mask <- mask(altCAN, qc) ``` --- # Crop and mask rasters ```r plot(alt_qc) ``` <img src="index_files/figure-html/unnamed-chunk-80-1.png" style="display: block; margin: auto;" /> --- # Set new raster projection We can use `projectRaster()` to transform the CRS of one spatial object to match another spatial object. ```r st_crs(sf_ca_cities_nad83) ``` ``` ## Coordinate Reference System: ## EPSG: 3347 ## proj4string: "+proj=lcc +lat_1=49 +lat_2=77 +lat_0=63.390675 +lon_0=-91.86666666666666 +x_0=6200000 +y_0=3000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs" ``` ```r alt_qc_nad83 <- projectRaster(alt_qc, 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 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs") ``` --- # Set new raster projection We can use `projectRaster()` to transform the CRS of one spatial object to match another spatial object. ```r plot(alt_qc_nad83) ``` <img src="index_files/figure-html/unnamed-chunk-82-1.png" style="display: block; margin: auto;" /> --- class: clear, middle #
Practice 1. Read with `sf` the [canadian cities CSV file](https://simplemaps.com/static/data/country-cities/ca/ca.csv). 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. Crop the elevation raster with those buffers (`sf::st_crop()`) and compute the mean per buffer --- # Stacking rasters Two possibilities: over variables or over time period. <img src="./assets/img/raster_pic.png" width="90%" style="display: block; margin: auto;" /> ??? Bricks vs Stacks --- # Stacking rasters Example from scratch ```r library(raster) # Create two empty rasters rs_var1 <- raster(ncol=10, nrow=10) rs_var2 <- raster(ncol=10, nrow=10) # Fill both with random values rs_var1[] <- runif(100) rs_var2[] <- runif(100) # Stack them st_vars <- stack(rs_var1, rs_var2) st_vars ``` ``` ## class : RasterStack ## dimensions : 10, 10, 100, 2 (nrow, ncol, ncell, nlayers) ## resolution : 36, 18 (x, y) ## extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax) ## crs : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 ## names : layer.1, layer.2 ## min values : 0.005568851, 0.006592623 ## max values : 0.9932329, 0.9891654 ``` --- # Coerce rasterstack or raster to dataframe ```r as.data.frame(st_vars, xy = TRUE) ``` ``` ## x y layer.1 layer.2 ## 1 -162 81 0.300138393 0.347314287 ## 2 -126 81 0.314687204 0.670251390 ## 3 -90 81 0.527173547 0.390061589 ## 4 -54 81 0.173962424 0.521274575 ## 5 -18 81 0.843360040 0.897462344 ## 6 18 81 0.040611388 0.710755422 ## 7 54 81 0.822329881 0.305853544 ## 8 90 81 0.108922448 0.308770588 ## 9 126 81 0.749261227 0.321085830 ## 10 162 81 0.752019061 0.249497145 ## 11 -162 63 0.233203975 0.417025461 ## 12 -126 63 0.773571847 0.979907134 ## 13 -90 63 0.443967437 0.511184320 ## 14 -54 63 0.437517612 0.172439877 ## 15 -18 63 0.615307346 0.035050276 ## 16 18 63 0.861475644 0.122374752 ## 17 54 63 0.575491242 0.749783191 ## 18 90 63 0.732974355 0.390325907 ## 19 126 63 0.166181299 0.087347549 ## 20 162 63 0.763771369 0.107741837 ## 21 -162 45 0.005568851 0.413725639 ## 22 -126 45 0.147631873 0.477088107 ## 23 -90 45 0.783283943 0.905040349 ## 24 -54 45 0.602983333 0.309525587 ## 25 -18 45 0.783841859 0.241424435 ## 26 18 45 0.653489037 0.183681450 ## 27 54 45 0.066448686 0.345308181 ## 28 90 45 0.370190436 0.707793848 ## 29 126 45 0.483124000 0.565369271 ## 30 162 45 0.226502897 0.699981068 ## 31 -162 27 0.934103163 0.880770787 ## 32 -126 27 0.232147107 0.397079693 ## 33 -90 27 0.426938216 0.624805124 ## 34 -54 27 0.147960207 0.487236475 ## 35 -18 27 0.357451751 0.656240061 ## 36 18 27 0.962355986 0.769876024 ## 37 54 27 0.421498507 0.436238035 ## 38 90 27 0.617186978 0.363403855 ## 39 126 27 0.188791232 0.619617144 ## 40 162 27 0.731675929 0.085550138 ## 41 -162 9 0.285505962 0.983396456 ## 42 -126 9 0.689455539 0.532673673 ## 43 -90 9 0.270820510 0.800885746 ## 44 -54 9 0.622043080 0.211926571 ## 45 -18 9 0.030688623 0.020714986 ## 46 18 9 0.457037790 0.550632767 ## 47 54 9 0.325140742 0.947222339 ## 48 90 9 0.747626565 0.364587458 ## 49 126 9 0.911774238 0.012390116 ## 50 162 9 0.839302338 0.592834600 ## 51 -162 -9 0.467755808 0.032522116 ## 52 -126 -9 0.238285254 0.115613834 ## 53 -90 -9 0.737491569 0.682802156 ## 54 -54 -9 0.675771538 0.705467062 ## 55 -18 -9 0.175756536 0.896276645 ## 56 18 -9 0.234526771 0.534642871 ## 57 54 -9 0.210108884 0.228574818 ## 58 90 -9 0.790257906 0.078019828 ## 59 126 -9 0.618018565 0.882332584 ## 60 162 -9 0.501333668 0.614874375 ## 61 -162 -27 0.033157347 0.510991049 ## 62 -126 -27 0.923703563 0.632586530 ## 63 -90 -27 0.401949014 0.140513252 ## 64 -54 -27 0.797342834 0.039361906 ## 65 -18 -27 0.947103777 0.017012870 ## 66 18 -27 0.058981298 0.808409599 ## 67 54 -27 0.746261590 0.463341615 ## 68 90 -27 0.554813607 0.354121953 ## 69 126 -27 0.073060868 0.938514855 ## 70 162 -27 0.935828489 0.684283321 ## 71 -162 -45 0.281746867 0.641507648 ## 72 -126 -45 0.093245816 0.913042879 ## 73 -90 -45 0.749709571 0.508319371 ## 74 -54 -45 0.260365228 0.117844829 ## 75 -18 -45 0.697727602 0.968789996 ## 76 18 -45 0.210735367 0.430471796 ## 77 54 -45 0.920472169 0.611271959 ## 78 90 -45 0.121831249 0.656451029 ## 79 126 -45 0.593781881 0.911441379 ## 80 162 -45 0.830561861 0.006592623 ## 81 -162 -63 0.860658494 0.356298496 ## 82 -126 -63 0.788011649 0.133621647 ## 83 -90 -63 0.662313491 0.077091290 ## 84 -54 -63 0.993232852 0.159367157 ## 85 -18 -63 0.318595531 0.152777357 ## 86 18 -63 0.626778644 0.283329941 ## 87 54 -63 0.384003210 0.484702571 ## 88 90 -63 0.708498249 0.989165384 ## 89 126 -63 0.180604049 0.616241620 ## 90 162 -63 0.864245176 0.450599653 ## 91 -162 -81 0.941628309 0.193843063 ## 92 -126 -81 0.370116698 0.683334720 ## 93 -90 -81 0.032062873 0.392572330 ## 94 -54 -81 0.497134914 0.507555191 ## 95 -18 -81 0.990798303 0.521776930 ## 96 18 -81 0.342713256 0.747875770 ## 97 54 -81 0.383560135 0.592604058 ## 98 90 -81 0.587500575 0.513813187 ## 99 126 -81 0.275187365 0.341539436 ## 100 162 -81 0.953163290 0.652676052 ``` ```r as.data.frame(rs_var1, xy = TRUE) ``` --- # Basic arithmetics ### Some arithmetic operations - Classic operators: `+`, `*`, `-` etc. - `rs_vars1 + rs_var2` - sum: `sum(st_vars)` - mean: `mean(st_vars)` - variance: `var(rs_var1)` - covariance: `cov(rs_var1, rs_var2)` - Frequecy values: `hist(rs_var1)` --- # Great online resources 1/3 ### Good tutorials for spatial data in R - [Raster analysis in R](https://mgimond.github.io/megug2017/) - [Geocomputation with R](https://geocompr.robinlovelace.net/intro.html) - [Spatial data in R](https://github.com/Pakillo/R-GIS-tutorial/blob/master/R-GIS_tutorial.md) - [Document par Nicolas Casajus (fr)](https://qcbs.ca/wiki/_media/gisonr.pdf) - http://r-spatial.org/ - [Tutorial on datacamp](https://www.datacamp.com/courses/spatial-analysis-in-r-with-sf-and-raster) - [R in space - Insileco](https://insileco.github.io/2018/04/14/r-in-space---a-series/) ### `sf` manipulations - [sf vignette #4](https://cran.r-project.org/web/packages/sf/vignettes/sf4.html) - [Geocomputation with R](https://geocompr.robinlovelace.net/attr.html) - [Attribute manipulations](https://insileco.github.io/2018/04/09/r-in-space---attribute-manipulations/) - [Tidy spatial data in R](http://strimas.com/r/tidy-sf/) --- # Great online resources 2/3 ### Maps in R - [Introduction to visualising spatial data in R](https://cran.r-project.org/doc/contrib/intro-spatial-rl.pdf) - [Geocomputation with R](https://geocompr.robinlovelace.net/adv-map.html) - [choropleth](https://cengel.github.io/rspatial/4_Mapping.nb.html) - [leaflet](https://rstudio.github.io/leaflet/) - [Mapview](https://r-spatial.github.io/mapview/index.html) - [tmap](https://cran.r-project.org/web/packages/tmap/vignettes/tmap-nutshell.html) - [plotly](https://plot.ly/python/maps/) - [Animated maps](https://insileco.github.io/2017/07/05/animations-in-r-time-series-of-erythemal-irradiance-in-the-st.-lawrence/) --- # Great online resources 3/3 ### Get free data - [free data at country level](http://www.diva-gis.org/gdata) - [Quebec free data](http://mffp.gouv.qc.ca/le-ministere/acces-aux-donnees-gratuites/) - [sdmpredictors](https://cran.r-project.org/web/packages/sdmpredictors/index.html) access to many datasets in R - [find more spatial data](https://freegisdata.rtwilson.com/) - [create shapefile on line](http://geojson.io/) - EPSG: [link1](http://spatialreference.org/); [link2](http://epsg.io/)