06 November 2018

In a nutshell

Packages

Resources

Why?

Mapping in Ecology?

  1. show where your plots are

  2. show how variables are distributed spatially

  3. spatial analyses

  4. many, many tools

R as a GIS?

  1. Workflow

  2. Quite efficient

  3. Interface

Workflow

R for everything (almost)

  1. import your data
  2. format your data
  3. analyze your data
  4. visualize your data
  5. export your data
  6. create your own functions/packages

Efficiency

  1. well-defined classes
  2. many formats (+ convert different file)
  3. operations on geometries
  4. statistical analyses
  5. visualization

Interface

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

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
  1. writeOGR / writeGDAL: to write spatial objects
  2. 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

Simple manipulations

plot(mysf[1])

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

Simple manipulations

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

Useful links to get data

Let’s practice

Installation

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)