title

author

date

output

NumeriLab; R maps and GIS

Arthur de Grandpré

25 janvier 2019

html_document

keep_md
true

Context

Most would agree that data science should always start and end with data visualization. It allows the scientist to understand the data, its structure and distribution, but also to connect with the public, by creating easily understandable content to facilitate communication.

In the previous workshops, we have learned how to:

Efficiently visualize data and results
Manipulate Data to facilitate its analysis
Generate data efficiently in order to replicate real life trends

... and all of this makes it hard to leave the R environment! Still, most of us still are leaving the comfort of R whenever we are faced with spatial data, splitting our workflow and creating different codes, different processes and ultimately resulting in a loss of time and efficiency.

While R still can't fully replace a full GIS suit, it is able to perform most tasks that are required on an almost daily basis for people working with spatial data.

R as a GIS and Mapping Software

To address this problem, multiple libraries have been made available in R in order to deal with spatial data, its visualisation, its manipulation and its analysis. the most notable are probably:

GIS
sp : allow the creation and manipulation of spatial objects (spatialpointsdataframes) + apply CRS and stransformations rgdal : An interface to access gdal, a spatial data processing library
rgeos : An interface to access geos, a spatial vector data processing library
GeoJSON : package to read GeoJSON data (single file container for spatial data)

mapping
base R : plot() can still do most of the job
ggplot2 : still able to do most of the job, prettier
ggmap : extends ggplot2 mapping functionalities and enables to retrieve google style static maps
ggsn : adds scale bars and north arrows to ggmap leaflet : allows to build interactive maps, google map style
ggsn : for on map scale and north arrow, but might not work as expected

other useful packages
tmap ; raster ; spdep ; cartography ; maptools ; gstat

Goals of the workshop

The main goal of this workshop is to break the mindset in which R is unable to deal with spatial data or spatial analysis.
It is important to keep in mind that this is only an introduction and that it is possible to do almost anything on R, and in many different ways.
Also, it might be even more important to keep in mind that in many situations, while R could do it, you might want to consider other solutions. Using R for everything is fun and can be practical, but sometimes it's just more tedious than necessary.

By the end of this workshop, I wish that frequent GIS tasks do not feel intimidating for those who attended, such as:

reading and visualising point data
transforming the CRS from data
loading a basemap from the internet
reading and visualising a vector file
reading and visualising a raster file
making simple maps
clipping a raster with a polygon
run an unsupervised classification
compare vegetation indexes (NDVI)
obtain zonal statistics on buffers
obtain landscape metrics

Exercice:

note: all data for this exercise is located in the repository under the "Data" sub-directories

#Step 1. Prepare R environment and load the libraries

rm(list=ls())

library(sp) # for functions coordinates, proj4string, CRS
library(rgdal) # for functions readOGR
library(ggmap)
library(ggsn) #scalebar and north arrow for ggplot and ggmap
library(leaflet)
library(mapview)#; webshot::install_phantomjs() has to be run only once

library(raster) # for functions raster, plotRGB#library(rgeos)
library(cluster)
library(SDMTools)
library(rgeos)

Step 2. reading and visualising point data + adding a Google type basemap

Multiples solutions are available to read point data in the R environment depending on your database format. The most common data format is text, such as .csv files. When working with .csv files, you will often have two columns referring to the coordinates of your data, which can be used create a spatial object with the package sp.

If you have a shapefile of your data, it will be imported directly as a spatial object by using the function readOGR() from the rgdal package

df = read.csv("Data/csv/data_picom_HT.csv",sep=";")
# this text file contains x and y coordinates and a random associated value
class(df)

## [1] "data.frame"

head(df)

##   zone parcelle             zone_eco      lat      long sup.zone..ha.
## 1   18      18A Jeune pinŠde blanche 46.34558 -72.58600         0.213
## 2    4       4A Jeune pinŠde blanche 46.34992 -72.57578         0.552
## 3   13      13A                mixte 46.34743 -72.58229         1.441
## 4   17      17A                mixte 46.34656 -72.58300         2.789
## 5   17      17B                mixte 46.34653 -72.58397         2.782
## 6   17      17C                mixte 46.34597 -72.58511         2.782
##   sup.parcelle..ha. tot_drymass aerial_drymass root_drymass n_trees
## 1              0.04        8861           6877         1984      95
## 2              0.04       11379           8896         2483      50
## 3              0.04        6035           4692         1343      51
## 4              0.04       13721          10766         2955      57
## 5              0.04        9865           7726         2140      40
## 6              0.04       10449           8200         2250      39
##   total_C n_living_trees mean_dhp max_dhp mean_height  dom_sp est_age  X
## 1    4735             84     14.9    29.0        14.1     PIB      45 NA
## 2    6236             50     23.5    33.9        13.6     PIB      45 NA
## 3    3146             38     16.6    28.9        13.6     BOP      45 NA
## 4    7651             42     22.2    79.1        16.9 PIB/PIG      45 NA
## 5    5397             29     23.5    42.4        18.3     PIB      45 NA
## 6    5727             34     22.8    49.2        15.1     PIB      45 NA

Help sheet for CRS selection in rgdal and sp @ https://www.nceas.ucsb.edu/~frazier/RSpatialGuides/OverviewCoordinateReferenceSystems.pdf

coordinates(df) = c("long","lat")
class(df) # at this point R recognizes the data as a spatial object

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

proj4string(df) = CRS("+init=epsg:4269") #this allows to specify the data CRS is NAD83

At this point, we have the same output as if we had used the function readOGR from rgdal to read a shapefile of the same information

summary(df)

## Object of class SpatialPointsDataFrame
## Coordinates:
##            min       max
## long -72.58601 -72.57278
## lat   46.34263  46.35089
## Is projected: FALSE 
## proj4string :
## [+init=epsg:4269 +proj=longlat +datum=NAD83 +no_defs +ellps=GRS80
## +towgs84=0,0,0]
## Number of points: 36
## Data attributes:
##       zone          parcelle                    zone_eco  sup.zone..ha.   
##  Min.   : 1.00   11A    : 1   PinŠde grise          :17   Min.   :0.1510  
##  1st Qu.: 7.00   12A    : 1   mixte                 : 5   1st Qu.:0.4895  
##  Median :17.00   12B    : 1   Vieille pinŠde blanche: 4   Median :1.5090  
##  Mean   :14.56   13A    : 1   Mixte surann‚         : 3   Mean   :1.3419  
##  3rd Qu.:20.00   14A    : 1   Jeune pinŠde blanche  : 2   3rd Qu.:1.8460  
##  Max.   :26.00   16A    : 1   PinŠde blanche_grise  : 2   Max.   :2.7890  
##                  (Other):30   (Other)               : 3                   
##  sup.parcelle..ha.  tot_drymass    aerial_drymass   root_drymass 
##  Min.   :0.04      Min.   : 4048   Min.   : 3159   Min.   : 889  
##  1st Qu.:0.04      1st Qu.: 6286   1st Qu.: 4902   1st Qu.:1384  
##  Median :0.04      Median : 7590   Median : 5922   Median :1668  
##  Mean   :0.04      Mean   : 8850   Mean   : 6927   Mean   :1924  
##  3rd Qu.:0.04      3rd Qu.:11498   3rd Qu.: 8996   3rd Qu.:2502  
##  Max.   :0.04      Max.   :20148   Max.   :15939   Max.   :4209  
##                                                                  
##     n_trees         total_C      n_living_trees    mean_dhp    
##  Min.   :18.00   Min.   : 2112   Min.   :17.0   Min.   :14.90  
##  1st Qu.:33.75   1st Qu.: 3306   1st Qu.:27.0   1st Qu.:18.40  
##  Median :40.00   Median : 4043   Median :33.5   Median :20.50  
##  Mean   :41.39   Mean   : 4782   Mean   :35.5   Mean   :21.43  
##  3rd Qu.:49.00   3rd Qu.: 6254   3rd Qu.:42.0   3rd Qu.:23.50  
##  Max.   :95.00   Max.   :11114   Max.   :84.0   Max.   :39.90  
##                                                                
##     max_dhp       mean_height        dom_sp      est_age     
##  Min.   :23.80   Min.   :12.80   PIG    :15   Min.   :45.00  
##  1st Qu.:29.80   1st Qu.:14.10   PIB    : 7   1st Qu.:45.00  
##  Median :37.65   Median :15.50   PIG/PIB: 3   Median :45.00  
##  Mean   :39.71   Mean   :16.23   Mixte  : 2   Mean   :51.94  
##  3rd Qu.:43.08   3rd Qu.:17.12   PIB/PIR: 2   3rd Qu.:55.00  
##  Max.   :79.10   Max.   :25.90   BOP    : 1   Max.   :70.00  
##                                  (Other): 6                  
##     X          
##  Mode:logical  
##  NA's:36       
##                
##                
##                
##                
##

It is very easy to transform the CRS of an object by using the spTransform function. it can be done my manually entering the datum, ellipsoid and projection, but the simplest way is to just feed the epsg code to the function. epsg codes are unique identifiers for CRS and can be bound easily on internet.

df <- spTransform(df, CRS("+init=epsg:4326"))  # transform to wgs84
summary(df)

## Object of class SpatialPointsDataFrame
## Coordinates:
##            min       max
## long -72.58601 -72.57278
## lat   46.34263  46.35089
## Is projected: FALSE 
## proj4string :
## [+init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84
## +towgs84=0,0,0]
## Number of points: 36
## Data attributes:
##       zone          parcelle                    zone_eco  sup.zone..ha.   
##  Min.   : 1.00   11A    : 1   PinŠde grise          :17   Min.   :0.1510  
##  1st Qu.: 7.00   12A    : 1   mixte                 : 5   1st Qu.:0.4895  
##  Median :17.00   12B    : 1   Vieille pinŠde blanche: 4   Median :1.5090  
##  Mean   :14.56   13A    : 1   Mixte surann‚         : 3   Mean   :1.3419  
##  3rd Qu.:20.00   14A    : 1   Jeune pinŠde blanche  : 2   3rd Qu.:1.8460  
##  Max.   :26.00   16A    : 1   PinŠde blanche_grise  : 2   Max.   :2.7890  
##                  (Other):30   (Other)               : 3                   
##  sup.parcelle..ha.  tot_drymass    aerial_drymass   root_drymass 
##  Min.   :0.04      Min.   : 4048   Min.   : 3159   Min.   : 889  
##  1st Qu.:0.04      1st Qu.: 6286   1st Qu.: 4902   1st Qu.:1384  
##  Median :0.04      Median : 7590   Median : 5922   Median :1668  
##  Mean   :0.04      Mean   : 8850   Mean   : 6927   Mean   :1924  
##  3rd Qu.:0.04      3rd Qu.:11498   3rd Qu.: 8996   3rd Qu.:2502  
##  Max.   :0.04      Max.   :20148   Max.   :15939   Max.   :4209  
##                                                                  
##     n_trees         total_C      n_living_trees    mean_dhp    
##  Min.   :18.00   Min.   : 2112   Min.   :17.0   Min.   :14.90  
##  1st Qu.:33.75   1st Qu.: 3306   1st Qu.:27.0   1st Qu.:18.40  
##  Median :40.00   Median : 4043   Median :33.5   Median :20.50  
##  Mean   :41.39   Mean   : 4782   Mean   :35.5   Mean   :21.43  
##  3rd Qu.:49.00   3rd Qu.: 6254   3rd Qu.:42.0   3rd Qu.:23.50  
##  Max.   :95.00   Max.   :11114   Max.   :84.0   Max.   :39.90  
##                                                                
##     max_dhp       mean_height        dom_sp      est_age     
##  Min.   :23.80   Min.   :12.80   PIG    :15   Min.   :45.00  
##  1st Qu.:29.80   1st Qu.:14.10   PIB    : 7   1st Qu.:45.00  
##  Median :37.65   Median :15.50   PIG/PIB: 3   Median :45.00  
##  Mean   :39.71   Mean   :16.23   Mixte  : 2   Mean   :51.94  
##  3rd Qu.:43.08   3rd Qu.:17.12   PIB/PIR: 2   3rd Qu.:55.00  
##  Max.   :79.10   Max.   :25.90   BOP    : 1   Max.   :70.00  
##                                  (Other): 6                  
##     X          
##  Mode:logical  
##  NA's:36       
##                
##                
##                
##                
##

at this point, the dataframe can be plotted in base R for visualisation...

plot(df)

..and made into a simple plot for visualisation.

plot(df,
     main="data points",
     axes=T,
     pch=21,
     bg=rev(heat.colors(5))[cut(df$tot_drymass,breaks=c(0,5000,7500,10000,12500,25000), labels = F)])

and the same can be done using ggplot2:

ggplot()+
  geom_point(data=as.data.frame(df), aes(long,lat, fill=tot_drymass), shape=21, size=2)+
  scale_fill_gradientn(colours=rev(heat.colors(5)))+
  ggtitle("data points")

So we can now see that we have spatial points with some values, but that doesn't give us much of an idea of what they are and their context. For adding context, a basemap is really important.

Note: Until recently, get_map from the package ggmap allowed to easily obtain basemaps from google and OpenStreetMaps. Sadly, google now requires an API key and billing information to access the basemaps, and it also broke the access to open access maps such as stamen. That makes loading a static basemap a bit harder. While it is still possible to access google maps and OpenStreetMaps, it is required to register an API key (generally free to obtain, but requires registration)

bbox(df)

##            min       max
## long -72.58601 -72.57278
## lat   46.34263  46.35089

map = get_stamenmap(bbox=c(left    = -72.6,
                           bottom  =  46.335,
                           right   = -72.565,
                           top     =  46.355),
                    zoom=14, maptype= "toner-lite")

# other basemap types from stamen:
#“terrain”, “terrain-background”, “terrain-labels”, “terrain-lines”, “toner”, “toner-2010”, “toner-2011”, “toner-background”, “toner-hybrid”, “toner-labels”, “toner-lines”, “toner-lite”, “watercolor”

Luckily, leaflet still does this very good for interactive maps.
By default, leaflet uses OpenStreetMap

# using the default values, leaflet uses OpenStreetMap
m1 <- leaflet(df) %>%
  addTiles()
m1

But it can be set to call other data sources, including google if an API key is registered.

m2 <- leaflet(df) %>%
  addProviderTiles("Stamen.Toner")
m2

It possesses different types of markers, such as pins.

m1 %>%
  addMarkers()

Or circles.

m1 %>%
  addCircleMarkers()

And those can be customized to display data values with continuous scales

pal = colorNumeric(
  palette = "RdYlGn",
  domain = df$tot_drymass
)

m1 %>%
  addCircleMarkers(color = "black", opacity=1, fillColor=~pal(tot_drymass), fillOpacity = 0.5, radius=8) %>%
  addLegend("bottomright",
            pal=pal,
            values=~tot_drymass,
            title= "drymass (kg / 0.04 ha")

Or by quantiles

qpal = colorQuantile("RdYlGn", df$tot_drymass, n = 4)

m1 %>%
  addCircleMarkers(color = "black", opacity=1, fillColor=~qpal(tot_drymass), fillOpacity = 0.5, radius=8) %>%
  addLegend("bottomright",
            pal=qpal,
            values=~tot_drymass,
            title= "Drymass Quantile") %>%
  addScaleBar(position="bottomleft")

And their size can also be fixed to represent the size of a variable

m1 %>%
  addCircles(color = "black", opacity=0.5, fillColor=~pal(tot_drymass), fillOpacity = 0.8, radius=df$n_trees) %>%
  addLegend("bottomright",
            pal=qpal,
            values=~tot_drymass,
            title= "Drymass Quantile") %>%
  addScaleBar(position="bottomleft")

And then this map can be saved to a png or jpg format using the mapview package

m1 = m1 %>%
  addCircles(stroke=F, opacity=1, fillColor=~pal(df$tot_drymass), fillOpacity = 0.8, radius=df$n_trees) %>%
  addLegend("bottomright",
            pal=qpal,
            values=~df$tot_drymass,
            title= "values") %>%
  addScaleBar()

mapshot(m1, file= "Figures/m1.png")

#Step 2. Loading and plotting vectors and rasters

often, we have polygons delimiting zones of interest for our data. they can be used for visualisation, but also in data analysis.

##adding polygon
campus= readOGR("Data/shapefiles/campus polygon.shp")

## OGR data source with driver: ESRI Shapefile 
## Source: "C:\Users\Arthur\Desktop\NumeriLab_RGIS\Data\shapefiles\campus polygon.shp", layer: "campus polygon"
## with 1 features
## It has 1 fields

campus= spTransform(campus, CRS("+init=epsg:4326"))
campus

## class       : SpatialPolygonsDataFrame 
## features    : 1 
## extent      : -72.58708, -72.57161, 46.34086, 46.35128  (xmin, xmax, ymin, ymax)
## coord. ref. : +init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## variables   : 1
## names       : Id 
## value       :  0

m1 = m1 %>% addPolygons(data=campus, stroke=T, fillOpacity = 0)
m1

When working with spatial data, it can often be useful to access remote sensing data to perform different types of analysis.

the rest of the workshop will be mostly about working with such data in R, based on a sample of Sentinel-2 imagery.

Then, let's load some data raster data that we can work with. To do so, the package "raster" contains most of the basic functions required for R to read and manage rasters. A raster is basically a spatially referenced matrix. where every value is associated with a spatial extent (resolution) and a coordinate taken from a georeference system (ex: WGS84).

# The raster function allows to read the most common raster formats, such as .TIFF
tr_r = raster("./Data/S2/True color.tiff")

# base R has some functions that allow for plotting of rasters
plot(tr_r, main="UQTR on 2018-09-13, Sentinel-2")

![](NumeriLab_RGIS_files/figure-html/loading and plotting rasters-1.png)

# something seems quite wrong with the colour, so let's inspect our raster
tr_r

## class       : RasterLayer 
## band        : 1  (of  3  bands)
## dimensions  : 212, 459, 97308  (nrow, ncol, ncell)
## resolution  : 10, 10  (x, y)
## extent      : -8081880, -8077290, 5834930, 5837050  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs 
## data source : C:\Users\Arthur\Desktop\NumeriLab_RGIS\Data\s2\True color.tiff 
## names       : True_color

# it seems we only have access to one of the 3 layers contained in tr_r. This is because the raster function is programmed to create a single layer object.
# still, we can look into better ways of diplaying this single band, such as grayscale
grayscale_colors = gray.colors(100,
                               start= 0.0,
                               end=1,
                               gamma=2.2,
                               alpha=NULL)
plot(tr_r, main="UQTR on 2018-09-13, red band grayscale, Sentinel-2", col=grayscale_colors)

# or load another band from the same source
tr_r2=raster("./Data/s2/True color.tiff", band=2)
tr_r2 # this shows we are now using band 2 of 3 instead of 1

## class       : RasterLayer 
## band        : 2  (of  3  bands)
## dimensions  : 212, 459, 97308  (nrow, ncol, ncell)
## resolution  : 10, 10  (x, y)
## extent      : -8081880, -8077290, 5834930, 5837050  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs 
## data source : C:\Users\Arthur\Desktop\NumeriLab_RGIS\Data\s2\True color.tiff 
## names       : True_color

# we need to use the brick function to create a RasterBrick, composed of multiple raster layers
tr_r = brick("./Data/s2/True color.tiff")
tr_r # we do have our 3 layers

## class       : RasterBrick 
## dimensions  : 212, 459, 97308, 3  (nrow, ncol, ncell, nlayers)
## resolution  : 10, 10  (x, y)
## extent      : -8081880, -8077290, 5834930, 5837050  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs 
## data source : C:\Users\Arthur\Desktop\NumeriLab_RGIS\Data\s2\True color.tiff 
## names       : True_color.1, True_color.2, True_color.3

tr_r=addLayer(tr_r, "Data/s2/B08.tiff")

tr_r=projectRaster(tr_r, crs=("+init=epsg:4326"))
plot(tr_r, col=grayscale_colors, main="UQTR on 2018-09-13, grayscale, Sentinel=2") # but the plot function does not plot them together. To do so, we need to specify each color bands, or use the function plotRGB

hist(tr_r)

plotRGB(tr_r, r=1,g=2,b=3,stretch="hist")

map1=viewRGB(tr_r, 1,2,3, map=m1, quantiles = c(0.05, 0.95))
map1

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

NumeriLab_RGIS.md

NumeriLab_RGIS.md

Context

R as a GIS and Mapping Software

Goals of the workshop

Exercice:

acknowledgements

additionnal ressources

Files

NumeriLab_RGIS.md

Latest commit

History

NumeriLab_RGIS.md

File metadata and controls

Context

R as a GIS and Mapping Software

Goals of the workshop

Exercice:

acknowledgements

additionnal ressources