Mapping in R

06 November 2018

In a nutshell

Packages

sf vector files
raster raster
mapview visualize

Resources

Why?

Mapping in Ecology?

show where your plots are
show how variables are distributed spatially
spatial analyses
many, many tools

R as a GIS?

Workflow
Quite efficient
Interface

Workflow

R for everything (almost)

import your data
format your data
analyze your data
visualize your data
export your data
create your own functions/packages

Efficiency

well-defined classes
many formats (+ convert different file)
operations on geometries
statistical analyses
visualization

Interface

C/C++: Geos, proj4, Gdal
JavaScript: Leaflet
otherGIS: RQGIS, rgrass7, RgoogleMaps
retrieve data: osmdata

How? An overview

Turning R into a powerful GIS

1- Define classes

sp
raster
sf

2- Topology operations

rgeos
sf

3- Import / export

rgdal
sf

Turning R into a powerful GIS

4- Analyses (few examples)

spdep
spatial
gstat
dismo

5- Creating and editing maps

graphics
ggplot
interactive plot / web-based plot
a blog post about this

Turning R into a powerful GIS

6- But…

Geo-referencing
Visualizing large spatial objects
Watershed analysis

Package `sp`

Classes

Classes / Functions	Contents
Points	list of points (set of coordinates)
SpatialPoints	list of points + CRS
SpatialPointsDataPoints	list of points + CRS + attribute table
Line	a line (set of coordinates)
Lines	list of lines
SpatialLines	list of lines + CRS
SpatialLinesDataFrame	list of lines + CRS + attribute table

Example: `SpatialPointsDataPoints`

long	lat	var1	var2
-80.95362	42.17096	-0.8852481	3.5288935
-80.36508	43.10166	1.7148959	6.6046439
-80.65781	42.39208	-0.2219311	0.2493644
-81.42152	42.09180	-2.5043391	6.3653082
-80.78941	43.26784	0.1766467	3.7825111
-81.10533	42.11308	-0.8977658	2.7027611

Example: `SpatialPointsDataFrame`

library(sp)

## Loading required package: methods

## 
## Attaching package: 'sp'

## The following object is masked from 'package:graphicsutils':
## 
##     compassRose

mysp <- SpatialPointsDataFrame(
  coords = mydata[,1:2],
  data = mydata[,3:4],
  proj4string = CRS(
    "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0"
  )
)

Example: `SpatialPointsDataFrame`

class(mysp)

## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"

isS4(mysp)

## [1] TRUE

slotNames(mysp)

## [1] "data"        "coords.nrs"  "coords"      "bbox"        "proj4string"

Example: `SpatialPointsDataFrame`

mysp@proj4string

## CRS arguments:
##  +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0

head(mysp@data)

##         var1      var2
## 1 -0.8852481 3.5288935
## 2  1.7148959 6.6046439
## 3 -0.2219311 0.2493644
## 4 -2.5043391 6.3653082
## 5  0.1766467 3.7825111
## 6 -0.8977658 2.7027611

Change projection: `spTransform`

(mysp2 <- spTransform(mysp, CRS=CRS("+proj=merc +ellps=GRS80")))

##            coordinates         var1      var2
## 1  (-9011715, 5157929) -0.885248073 3.5288935
## 2  (-8946200, 5298253)  1.714895877 6.6046439
## 3  (-8978786, 5191078) -0.221931147 0.2493644
## 4  (-9063802, 5146090) -2.504339144 6.3653082
## 5  (-8993436, 5323533)  0.176646727 3.7825111
## 6  (-9028604, 5149271) -0.897765848 2.7027611
## 7  (-8985002, 5428508) -1.004612170 6.5370860
## 8  (-8987188, 5403104)  0.081480485 5.6064616
## 9  (-9112177, 5406733)  0.887519582 2.6685859
## 10 (-9121702, 5410634)  0.175649243 4.2287652
## 11 (-8916394, 5286647) -0.143991329 9.3165423
## 12 (-8908251, 5423223)  0.014687077 1.1272343
## 13 (-9103152, 5331383) -0.448895976 3.3232480
## 14 (-9121535, 5338572)  0.561188239 0.2582339
## 15 (-9126676, 5419720) -0.825143542 2.2148579
## 16 (-9032011, 5147282)  0.993920452 1.3432549
## 17 (-9035231, 5326703) -0.003455062 7.2756258
## 18 (-9080838, 5387184) -0.660376744 7.6992890
## 19 (-8994340, 5360932)  0.276733107 8.9131325
## 20 (-9047506, 5136935)  0.339732272 9.9522539

Package `raster`

Classes

0- SpatialGrid / SpatialPixel (package `sp`)

3- RasterLayer (package `raster`)

4- RasterStack (package `raster`)

5- RasterBrick (package `raster`)

Example: `raster`

library(raster)

## 
## Attaching package: 'raster'

## The following object is masked from 'package:magrittr':
## 
##     extract

ras1 <- raster(matrix(runif(100*100,0,10),ncol=100,nrow=100),
    crs=CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"),
    xmn=-82, xmx=+80, ymn=+42, ymx=+44)

Example: `raster`

ras1

## class       : RasterLayer 
## dimensions  : 100, 100, 10000  (nrow, ncol, ncell)
## resolution  : 1.62, 0.02  (x, y)
## extent      : -82, 80, 42, 44  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : in memory
## names       : layer 
## values      : 0.001514384, 9.999197  (min, max)

Retrieving free GIS data: `getData`

head(getData("ISO3"))

##   ISO3                  NAME
## 1  AFG           Afghanistan
## 2  XAD Akrotiri and Dhekelia
## 3  ALA                 Åland
## 4  ALB               Albania
## 5  DZA               Algeria
## 6  ASM        American Samoa

Example: `getData`

## Country level:
mapBEL0 <- getData(name="GADM", country="BEL", path="./assets", level=0)
mapBEL1 <- getData(name="GADM", country="BEL", path="./assets", level=1)
tminW <- getData(name="worldclim", var="tmin", res=10, path="./assets")
mapBEL0

## class       : SpatialPolygonsDataFrame 
## features    : 1 
## extent      : 2.555356, 6.40787, 49.49722, 51.50382  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## variables   : 2
## names       : GID_0,  NAME_0 
## value       :   BEL, Belgium

Package `rgdal`

Drivers

library(rgdal)

## rgdal: version: 1.3-6, (SVN revision 773)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.3.2, released 2018/09/21
##  Path to GDAL shared files: /usr/share/gdal
##  GDAL binary built with GEOS: TRUE 
##  Loaded PROJ.4 runtime: Rel. 5.2.0, September 15th, 2018, [PJ_VERSION: 520]
##  Path to PROJ.4 shared files: (autodetected)
##  Linking to sp version: 1.3-1

writeOGR / writeGDAL: to write spatial objects
readOGR/readGDAL: to read spatial files

Drivers

head(ogrDrivers())

##         name                     long_name write  copy isVector
## 1 AeronavFAA                   Aeronav FAA FALSE FALSE     TRUE
## 2 AmigoCloud                    AmigoCloud  TRUE FALSE     TRUE
## 3     ARCGEN             Arc/Info Generate FALSE FALSE     TRUE
## 4     AVCBin      Arc/Info Binary Coverage FALSE FALSE     TRUE
## 5     AVCE00 Arc/Info E00 (ASCII) Coverage FALSE FALSE     TRUE
## 6        BNA                     Atlas BNA  TRUE FALSE     TRUE

Export

writeOGR(mysp, dsn="./assets", layer="mypoints",
    driver="ESRI Shapefile", overwrite_layer=TRUE)

Import

mysp2 <- readOGR(dsn="assets/", layer="mypoints")

## OGR data source with driver: ESRI Shapefile 
## Source: "/home/kevcaz/Github/Tutorials/mapsWithR/docs/assets", layer: "mypoints"
## with 20 features
## It has 2 fields

## Roads find on http://www.diva-gis.org/Data
canroads <- readOGR(dsn="assets/", layer="CAN_roads")

## OGR data source with driver: ESRI Shapefile 
## Source: "/home/kevcaz/Github/Tutorials/mapsWithR/docs/assets", layer: "CAN_roads"
## with 16915 features
## It has 5 fields

NB: there are function to import/export raster in raster

Package `rgeos`

Load the package

library(rgeos)

## rgeos version: 0.3-28, (SVN revision 572)
##  GEOS runtime version: 3.7.0-CAPI-1.11.0 673b9939 
##  Linking to sp version: 1.3-1 
##  Polygon checking: TRUE

Load the package

buf <- gBuffer(mapBEL0, width=0.5)

## Warning in gBuffer(mapBEL0, width = 0.5): Spatial object is not projected;
## GEOS expects planar coordinates

diff <- gDifference(gBuffer(mapBEL0, width=0.5), gBuffer(mapBEL0, width=0.1))

## Warning in gBuffer(mapBEL0, width = 0.5): Spatial object is not projected;
## GEOS expects planar coordinates

## Warning in gBuffer(mapBEL0, width = 0.1): Spatial object is not projected;
## GEOS expects planar coordinates

package `sf`

Simple features?

The package sf […] aims at succeeding sp in the long term.

Simple features or simple feature access refers to a formal standard (ISO 19125-1:2004) that describes how objects in the real world can be represented in computers, with emphasis on the spatial geometry of these objects

https://cran.r-project.org/web/packages/sf/vignettes/sf1.html

Importing spatial object

library(sf)

## Linking to GEOS 3.6.2, GDAL 2.3.2, proj.4 5.1.0

mysf <- st_read(dsn="assets/mypoints.shp")

## Reading layer `mypoints' from data source `/home/kevcaz/Github/Tutorials/mapsWithR/docs/assets/mypoints.shp' using driver `ESRI Shapefile'
## Simple feature collection with 20 features and 2 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: -81.98633 ymin: 42.03051 xmax: -80.02418 ymax: 43.95302
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs

Importing spatial object

print(mysf)

## Simple feature collection with 20 features and 2 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: -81.98633 ymin: 42.03051 xmax: -80.02418 ymax: 43.95302
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
## First 10 features:
##           var1      var2                   geometry
## 1  -0.88524807 3.5288935 POINT (-80.95362 42.17096)
## 2   1.71489588 6.6046439 POINT (-80.36508 43.10166)
## 3  -0.22193115 0.2493644 POINT (-80.65781 42.39208)
## 4  -2.50433914 6.3653082  POINT (-81.42152 42.0918)
## 5   0.17664673 3.7825111 POINT (-80.78941 43.26784)
## 6  -0.89776585 2.7027611 POINT (-81.10533 42.11308)
## 7  -1.00461217 6.5370860 POINT (-80.71365 43.95302)
## 8   0.08148049 5.6064616 POINT (-80.73328 43.78793)
## 9   0.88751958 2.6685859 POINT (-81.85608 43.81154)
## 10  0.17564924 4.2287652 POINT (-81.94164 43.83691)

Importing spatial object

mysf2 <- st_as_sf(mysp)

Simple manipulations

plot(mysf[1])

plot(st_buffer(mysf[1], 0.1))

Simple manipulations

## Warning in st_buffer.sfc(st_geometry(x), dist, nQuadSegs): st_buffer does
## not correctly buffer longitude/latitude data

## dist is assumed to be in decimal degrees (arc_degrees).

`sf` or `sp`?

Package `mapview`

Import package

library(mapview)

NB: it uses the leaflet package which uses sp.

1 line of code examples

mapview(mysp, cex = 'var2')@map

Quick examples

mapview(mapBEL1)@map

More (using `leaflet`)

Editing a map

A very basic map – Shapefile

plot(mapBEL0)

A very basic map – Raster

class(tminW)

## [1] "RasterStack"
## attr(,"package")
## [1] "raster"

A very basic map – Raster

plot(tminW)

Customize a map – Shapefile

plot(mapBEL0, border='grey15', col='#E6E6E6', lwd=1.6)
plot(mapBEL1, lty=2, lwd=0.9, add=T)
points(4.3513, 50.8471, pch=19, col="#27df9d")
text(4.3513, 50.8471, text="Brussel", pos=3)

Customize a map – Shapefile

Resources

Useful links to use R as a GIS

task view: https://cran.r-project.org/web/views/Spatial.html; also following to this https://cran.r-project.org/web/views/

install.packages("ctv")
install.views("Spatial")

Useful links to get data

Let’s practice

Installation

1- Have a look at https://github.com/r-spatial/sf

2- Install R packages

install.packages(c("sf", "sp", "rgeos", "rgdal", "raster", "mapview"))

Alexis’s data

tmp <- read.csv('assets/Feuil2.csv')
head(tmp)

##      plot      lat       lon     elev
## 1 M-FT-01 45.44728 -71.10890 749.4448
## 2 M-FB-01 45.44738 -71.11548 908.9915
## 3 M-FM-01 45.44701 -71.11367 848.7617
## 4 M-FB-02 45.45012 -71.11455 909.0892
## 5 M-FM-02 45.44985 -71.11202 833.5641
## 6 M-FT-02 45.44917 -71.10891 738.8082

alex <- SpatialPointsDataFrame(
  tmp[c('lon', 'lat')],
  data = tmp[c('plot', 'elev')],
  proj4string = CRS('+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0')
)
head(alex)

##      plot     elev
## 1 M-FT-01 749.4448
## 2 M-FB-01 908.9915
## 3 M-FM-01 848.7617
## 4 M-FB-02 909.0892
## 5 M-FM-02 833.5641
## 6 M-FT-02 738.8082

plot(alex)

mapview(alex)

coord <- coordinates(alex)
# altCAN2 <- raster::getData(name="SRTM", lon = coord[1,1], lat = coord[1,2], path="assets/")
altCAN <- raster::getData(name="alt", country = "CAN", path="assets/")
bouCAN2 <- raster::getData(country='CAN', level=2, path="assets/")
bouCAN2_Q <- bouCAN2[bouCAN2@data$NAME_1 == "Québec",]

par(mar=rep(0,4))
plot(bouCAN2_Q)

al_ov <- over(alex, bouCAN2_Q)
head(al_ov)

##   GID_0 NAME_0    GID_1 NAME_1 NL_NAME_1       GID_2    NAME_2 VARNAME_2
## 1   CAN Canada CAN.11_1 Québec      <NA> CAN.11.54_1 Le Granit      <NA>
## 2   CAN Canada CAN.11_1 Québec      <NA> CAN.11.54_1 Le Granit      <NA>
## 3   CAN Canada CAN.11_1 Québec      <NA> CAN.11.54_1 Le Granit      <NA>
## 4   CAN Canada CAN.11_1 Québec      <NA> CAN.11.54_1 Le Granit      <NA>
## 5   CAN Canada CAN.11_1 Québec      <NA> CAN.11.54_1 Le Granit      <NA>
## 6   CAN Canada CAN.11_1 Québec      <NA> CAN.11.54_1 Le Granit      <NA>
##   NL_NAME_2                       TYPE_2                    ENGTYPE_2 CC_2
## 1      <NA> Regional County Municipality Regional County Municipality   30
## 2      <NA> Regional County Municipality Regional County Municipality   30
## 3      <NA> Regional County Municipality Regional County Municipality   30
## 4      <NA> Regional County Municipality Regional County Municipality   30
## 5      <NA> Regional County Municipality Regional County Municipality   30
## 6      <NA> Regional County Municipality Regional County Municipality   30
##     HASC_2
## 1 CA.QC.GR
## 2 CA.QC.GR
## 3 CA.QC.GR
## 4 CA.QC.GR
## 5 CA.QC.GR
## 6 CA.QC.GR

id <- bouCAN2_Q@data$NAME_2 %in% al_ov$NAME_2
bouCAN2_Q@data$NAME_2[id]

## [1] "Le Granit"

par(mar=rep(0,4))
plot(bouCAN2_Q, lwd=.4)
plot(bouCAN2_Q[id,], add=T, border = NA, col = 2)

ra_elv <- rasterize(bouCAN2_Q[id,], crop(altCAN, bouCAN2_Q[id,]@bbox), mask=TRUE)
ra_elv

## class       : RasterLayer 
## dimensions  : 87, 120, 10440  (nrow, ncol, ncell)
## resolution  : 0.008333333, 0.008333333  (x, y)
## extent      : -71.39167, -70.39167, 45.23333, 45.95833  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
## data source : in memory
## names       : CAN_msk_alt 
## values      : 243, 1103  (min, max)

par(mar=rep(0,4))
plot(ra_elv)
contour(ra_elv, add=T)

buf <- gBuffer(alex, FALSE, width = .001)

## Warning in gBuffer(alex, FALSE, width = 0.001): Spatial object is not
## projected; GEOS expects planar coordinates

buf

## class       : SpatialPolygons 
## features    : 1 
## extent      : -71.11806, -71.1066, 45.44112, 45.46517  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0

plot(buf, col ="grey50")
plot(alex, add = TRUE, pch=20, col="grey20")

altCAN2 <- raster('assets/srtm_22_03.tif')
ra_elv2 <- rasterize(bouCAN2_Q[id,], crop(altCAN2, bouCAN2_Q[id,]@bbox), mask = TRUE)
alex@data$totalC <- runif(nrow(alex@data), 5, 50)
alex@data$categ <- rep(1:3, 10)
alex@data$elv2 <- extract(ra_elv, alex)
alex@data$elv3 <- extract(ra_elv2, alex)
head(alex@data)

##      plot     elev   totalC categ elv2 elv3
## 1 M-FT-01 749.4448 10.49777     1  802  723
## 2 M-FB-01 908.9915 11.29750     2  802  910
## 3 M-FM-01 848.7617 48.26493     3  802  834
## 4 M-FB-02 909.0892 18.43345     1  862  886
## 5 M-FM-02 833.5641 24.32812     2  802  815
## 6 M-FT-02 738.8082 37.34865     3  802  737

par(las=1)
palette(c('#c17abb', '#38e2d2', '#dadd65'))
plot(t(alex@bbox), type = 'n')
plot(alex, add=T, cex = log(alex@data$totalC), pch = 19, col = alex@data$categ)

In a nutshell

Packages

Resources

Why?

Mapping in Ecology?

R as a GIS?

Workflow

R for everything (almost)

Efficiency

Interface

How? An overview

Turning R into a powerful GIS

1- Define classes

2- Topology operations

3- Import / export

Turning R into a powerful GIS

4- Analyses (few examples)

5- Creating and editing maps

Turning R into a powerful GIS

6- But…

Package sp

Classes

Example: SpatialPointsDataPoints

Example: SpatialPointsDataFrame

Example: SpatialPointsDataFrame

Example: SpatialPointsDataFrame

Change projection: spTransform

Package raster

Classes

0- SpatialGrid / SpatialPixel (package sp)

3- RasterLayer (package raster)

4- RasterStack (package raster)

5- RasterBrick (package raster)

Example: raster

Example: raster

Retrieving free GIS data: getData

Example: getData

Package rgdal

Drivers

Drivers

Export

Import

Package rgeos

Load the package

Load the package

package sf

Simple features?

Importing spatial object

Importing spatial object

Importing spatial object

Simple manipulations

Simple manipulations

Simple manipulations

sf or sp?

Package mapview

Import package

1 line of code examples

Quick examples

More (using leaflet)

Editing a map

A very basic map – Shapefile

A very basic map – Raster

A very basic map – Raster

Customize a map – Shapefile

Customize a map – Shapefile

Resources

Useful links to use R as a GIS

Useful links to get data

Let’s practice

Installation

1- Have a look at https://github.com/r-spatial/sf

2- Install R packages

Alexis’s data

Package `sp`

Example: `SpatialPointsDataPoints`

Example: `SpatialPointsDataFrame`

Example: `SpatialPointsDataFrame`

Example: `SpatialPointsDataFrame`

Change projection: `spTransform`

Package `raster`

0- SpatialGrid / SpatialPixel (package `sp`)

3- RasterLayer (package `raster`)

4- RasterStack (package `raster`)

5- RasterBrick (package `raster`)

Example: `raster`

Example: `raster`

Retrieving free GIS data: `getData`

Example: `getData`

Package `rgdal`

Package `rgeos`

package `sf`

`sf` or `sp`?

Package `mapview`

More (using `leaflet`)