RGIS 3.0

Context

Most would agree that data science should always start and end with data visualization. It allows scientists to understand the data structure and distribution, but also to connect with the public, by creating easily understandable content that facilitates 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! Most of us still are leaving the comfort of R whenever we are faced with spatial data, splitting our workflows and creating different scripts, different processes and ultimately resulting in a loss of time and efficiency.

If you’re willing to dig deep enough, you’ll find out that R can be used as part of a high performance GIS workflow for many applications, especially when combined with some external programs such as QGIS, GRASS and GDAL.

The main goal of this workshop is to break the mindset in which R is bad at dealing 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 with 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.

Version 3.0

Previous versions of this workshop were mostly focused on understanding basic data structure for spatial objects from the “sp” package, and how they can be used for basic mapping and building interactive maps using leaflet.

This version will shift towards a more up-do-date approach to spatial data in R using the “sf” package instead of “sp”. sf is faster, integrates more functions, and is compatible with tidy writing. This workshop will also put more focus on producing actually publishable static maps, presenting interactive maps as a collaboration tool.

R as a GIS and mapping software

Multiple libraries have been made available in R in order to deal with spatial data, its visualization, its manipulation and its analysis.
The most notable are probably the following:

GIS
sp : allow the creation and manipulation of spatial objects (spatialpointsdataframes) + apply CRS and transformations sf : alternative environment to sp, integrating more functions (rgeos and rgdal) and compatible with tidy writing rgdal : An interface to access gdal, a spatial data processing library
rgeos : An interface to access geos, a spatial vector data processing library raster : Allows to read and manipulate raster data. Very important!

mapping
ggmap : extends ggplot2 mapping functionalities and is able to retrieve basemaps
leaflet : allows to build interactive maps, Google map style, easy to share
ggsn : add North arrow and scalebar to ggmap

Exercices

NOTE: All data for this exercise are located in the repository under the “Data” sub-directory, available here: https://github.com/arthurdegrandpre/NumeriLab_RGIS/tree/2021

Step 1. Prepare the R environment and install/load the libraries

install.packages("sp")
install.packages("sf")
install.packages("tidyverse")
install.packages("ggmap")
install.packages("leaflet")
install.packages("mapview")
install.packages("raster")
install.packages("leafsync")
install.packages("htmlwidgets")
install.packages("ggsn")
install.packages("rgdal", type = "source") # requires an installation of Rtools, see https://cran.r-project.org/bin/windows/Rtools/
webshot::install_phantomjs() # install only once
rm(list=ls())

library("sf")
library("tidyverse")
library("ggmap")
library("leaflet")
library("mapview")
library("raster")
library("leafsync")
library("htmlwidgets")
library("sp")
library("ggmap")
library("ggsn")
library("rgdal")

Step 2. Reading and visualising point data + adding a 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 packages sf or sp.

If you have a shapefile, it can be imported directly as a spatial object by using the function st_read from sf or readOGR() from the rgdal package. The later will generate an sp compatible object.

df = read.csv("../Data/csv/data_picom_HT.csv",sep=";") # contains forestry data for the UQTR campus, with lon/lat coordinates

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

All functions from sf are called with the prefix “st_”. “st_as_sf()” can be used to transform a foreign object into a sf object, followed by “st_set_crs()” to define coordinate system. The data is in NAD83, which corresponds to epsg code 4269.

dfs = df %>% 

  st_as_sf(coords = c("lon","lat")) %>% 
  st_set_crs("epsg:4269")

dfs
## Simple feature collection with 36 features and 16 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -72.58601 ymin: 46.34263 xmax: -72.57278 ymax: 46.35089
## Geodetic CRS:  NAD83
## First 10 features:
##    zone parcelle             zone_eco sup_zone_ha sup_parc_ha tot_drymass
## 1    18      18A jeune pinede blanche       0.213        0.04        8861
## 2     4       4A jeune pinede blanche       0.552        0.04       11379
## 3    13      13A                mixte       1.441        0.04        6035
## 4    17      17A                mixte       2.789        0.04       13721
## 5    17      17B                mixte       2.782        0.04        9865
## 6    17      17C                mixte       2.782        0.04       10449
## 7    17      17D                mixte       2.782        0.04        6449
## 8    20      20A        mixte suranne       1.761        0.04       13441
## 9    20      20B        mixte suranne       1.761        0.04       20148
## 10   20      20C        mixte suranne       1.761        0.04       13799
##    aerial_drymass root_drymass n_trees total_C n_living_trees mean_dhp max_dhp
## 1            6877         1984      95    4735             84     14.9    29.0
## 2            8896         2483      50    6236             50     23.5    33.9
## 3            4692         1343      51    3146             38     16.6    28.9
## 4           10766         2955      57    7651             42     22.2    79.1
## 5            7726         2140      40    5397             29     23.5    42.4
## 6            8200         2250      39    5727             34     22.8    49.2
## 7            5023         1426      49    3425             41     17.2    35.2
## 8           10583         2858      28    7375             27     23.8    55.4
## 9           15939         4209      18   11114             17     39.9    72.2
## 10          10865         2934      33    7316             32     22.6    68.8
##    mean_height  dom_sp est_age                   geometry
## 1         14.1     PIB      45   POINT (-72.586 46.34558)
## 2         13.6     PIB      45 POINT (-72.57578 46.34992)
## 3         13.6     BOP      45 POINT (-72.58229 46.34743)
## 4         16.9 PIB/PIG      45   POINT (-72.583 46.34656)
## 5         18.3     PIB      45 POINT (-72.58397 46.34653)
## 6         15.1     PIB      45 POINT (-72.58511 46.34597)
## 7         13.0 PIB/ERR      45  POINT (-72.58475 46.3465)
## 8         17.1   Mixte      50 POINT (-72.58601 46.34511)
## 9         20.0   Mixte      70 POINT (-72.58511 46.34481)
## 10        16.5 Feuillu      55  POINT (-72.5843 46.34298)

It is very easy to transform the CRS of an object by using the st_transform function. It can be done my manually by 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 found easily online. Let’s transform to WGS84 (EPSG 4326).

dfs = dfs %>% 
  st_transform("epsg:4326")

dfs
## Simple feature collection with 36 features and 16 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -72.58601 ymin: 46.34263 xmax: -72.57278 ymax: 46.35089
## Geodetic CRS:  WGS 84
## First 10 features:
##    zone parcelle             zone_eco sup_zone_ha sup_parc_ha tot_drymass
## 1    18      18A jeune pinede blanche       0.213        0.04        8861
## 2     4       4A jeune pinede blanche       0.552        0.04       11379
## 3    13      13A                mixte       1.441        0.04        6035
## 4    17      17A                mixte       2.789        0.04       13721
## 5    17      17B                mixte       2.782        0.04        9865
## 6    17      17C                mixte       2.782        0.04       10449
## 7    17      17D                mixte       2.782        0.04        6449
## 8    20      20A        mixte suranne       1.761        0.04       13441
## 9    20      20B        mixte suranne       1.761        0.04       20148
## 10   20      20C        mixte suranne       1.761        0.04       13799
##    aerial_drymass root_drymass n_trees total_C n_living_trees mean_dhp max_dhp
## 1            6877         1984      95    4735             84     14.9    29.0
## 2            8896         2483      50    6236             50     23.5    33.9
## 3            4692         1343      51    3146             38     16.6    28.9
## 4           10766         2955      57    7651             42     22.2    79.1
## 5            7726         2140      40    5397             29     23.5    42.4
## 6            8200         2250      39    5727             34     22.8    49.2
## 7            5023         1426      49    3425             41     17.2    35.2
## 8           10583         2858      28    7375             27     23.8    55.4
## 9           15939         4209      18   11114             17     39.9    72.2
## 10          10865         2934      33    7316             32     22.6    68.8
##    mean_height  dom_sp est_age                   geometry
## 1         14.1     PIB      45   POINT (-72.586 46.34558)
## 2         13.6     PIB      45 POINT (-72.57578 46.34992)
## 3         13.6     BOP      45 POINT (-72.58229 46.34743)
## 4         16.9 PIB/PIG      45   POINT (-72.583 46.34656)
## 5         18.3     PIB      45 POINT (-72.58397 46.34653)
## 6         15.1     PIB      45 POINT (-72.58511 46.34597)
## 7         13.0 PIB/ERR      45  POINT (-72.58475 46.3465)
## 8         17.1   Mixte      50 POINT (-72.58601 46.34511)
## 9         20.0   Mixte      70 POINT (-72.58511 46.34481)
## 10        16.5 Feuillu      55  POINT (-72.5843 46.34298)

At this point, R recognizes the data as spatial data and will plot it as such if prompted to. The plot can then be manipulated into a simple data representation using base R. (More interesting with polygons… see https://r-spatial.github.io/sf/articles/sf5.html)

plot(dfs)
## Warning: plotting the first 9 out of 16 attributes; use max.plot = 16 to plot
## all

plot(dfs["tot_drymass"],
     key.pos = 1, #(1=below, 2=left, 3=above and 4=right)
     axes = T,
     key.width = lcm(1.3),
     key.length = 1.0,
     ) 

plot(dfs["tot_drymass"],
     key.pos = 1, #(1=below, 2=left, 3=above and 4=right)
     axes = T,
     key.width = lcm(1.3),
     key.length = 1.0,
     breaks = "jenks")

and the same can be done using ggplot2:

ggplot(dfs)+
  geom_sf(aes(fill=tot_drymass), shape=21, size=3)+
  scale_fill_viridis_c()

We can now see that we have spatial points with some values, but that does not 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 (which can generate costs for very heavy users) .

https://www.r-bloggers.com/geocoding-with-ggmap-and-the-google-api/ for Google API

# download a basemap for our data using get_stamenmap, and the bounding box in "sp" format. adding a small buffer around the bounding box prevents the basemap from being too limited in it's spatial extent

map = get_stamenmap(bbox=bbox(as_Spatial(dfs))+c(-0.001,-0.001,0.001,0.001),
                  maptype = "terrain",
                  zoom = 15)

# 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”

# using ggmap with the ggplot2 syntax allows to create complete static maps with relative ease

map1 = ggmap(map) +
  geom_sf(data = dfs,
          aes(fill = tot_drymass),
          size = 3,
          shape = 21,
          inherit.aes=F
          ) +
  labs(
    title = "Estimated forest drymass on UQTR campus, 2016",
    # subtitle = "Subtitle", 
    caption = "Basemap = Stamen Terrain, CRS = EPSG:4326 WGS84, produced 2021-02-02",
    tag = "a)",
    fill = "Tot. Drymass (kg)",
    x = "Longitude",
    y = "Latitude"
    )+
  theme(
    axis.text.x = element_text(angle = 45, hjust=1),
    axis.text.y = element_text(angle = 45, hjust=1),
    plot.caption = element_text(hjust = 0.5),
    plot.title.position = "plot"
    )+
  north(dfs,scale=0.2)+
  scalebar(dfs, dist = 250, dist_unit = "m",
             transform = TRUE, model = "WGS84", st.bottom = TRUE, st.dist = 0.05)
  

map1

ggsave("../Figures/test01.jpg", plot=map1)

We are planning a more advanced workshop about R and static maps in the future, but the focus of this workshop is about interactive maps, more specifically using the Leaflet package.

Step 3. Using leaflet to create interactive maps

Luckily, leaflet is still able to call basemaps without registering an API key.
By default, leaflet uses OpenStreetMap.

# using the default values, leaflet uses OpenStreetMap
leaflet(dfs) %>%
  addTiles()

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

leaflet(dfs) %>%
  addProviderTiles("Stamen.Terrain")

It possesses different types of markers, such as pins.

leaflet(dfs) %>%
  addTiles() %>%
  addMarkers()

Or circles.

leaflet(dfs) %>%
  addTiles() %>%
  addCircleMarkers()

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

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

leaflet(dfs) %>%
  addTiles() %>%
  addCircleMarkers(color = "black", opacity=0.5, 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", dfs$tot_drymass, n = 4)

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

Marker size can also be fixed to represent the size of a variable.

leaflet(dfs) %>%
  addTiles() %>%
  addCircles(color = "black", opacity=0.5, fillColor=~pal(tot_drymass), fillOpacity = 0.8, radius=dfs$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 = leaflet(dfs) %>%
  addTiles() %>%
  addCircles(stroke=T,color="black", opacity=0.5, fillColor=~pal(dfs$tot_drymass), fillOpacity = 0.8, radius=dfs$n_trees) %>%
  addLegend("bottomright",
            pal=qpal,
            values=~dfs$tot_drymass,
            title= "Drymass quantile") %>%
  addScaleBar(position="bottomleft")

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

Step 4. Loading and plotting vectors and rasters

Often, we have polygons delimiting zones of interest for our data. They can be used for visualization, but also in data analysis.

##adding polygon
campus = st_read("../Data/shapefiles/campus polygon.shp") %>%
  st_transform("epsg:4326")
## Reading layer `campus polygon' from data source `C:\Users\Arthur\Documents\rscripts\NumeriLab_RGIS\Data\shapefiles\campus polygon.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 1 field
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -8080357 ymin: 5804215 xmax: -8078634 ymax: 5805890
## Projected CRS: World_Mercator
leaflet(dfs) %>%
  addTiles() %>%
  addCircles(stroke=T, color="black", opacity=0.5, fillColor=~pal(dfs$tot_drymass), fillOpacity = 0.8, radius=dfs$n_trees) %>%
  addLegend("bottomright",
            pal=qpal,
            values=~tot_drymass,
            title= "values") %>%
  addScaleBar(position="bottomleft") %>%
  addPolygons(data=campus, stroke=T, fillOpacity = 0)

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 raster 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 each value is associated with a spatial extent (resolution) and a coordinate taken from a geodetic 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")
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded ellps WGS 84 in CRS definition: +proj=merc +a=6378137
## +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
## +wktext +no_defs +type=crs
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum World Geodetic System 1984 in CRS definition
# base R has some functions that allow to plot rasters
plot(tr_r, main="UQTR on 2018-09-13, Sentinel-2")

# something seems quite wrong with the color, 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)
## crs        : +proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs
## source     : C:/Users/Arthur/Documents/rscripts/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 displaying 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)
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded ellps WGS 84 in CRS definition: +proj=merc +a=6378137
## +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
## +wktext +no_defs +type=crs
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum World Geodetic System 1984 in CRS definition
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)
## crs        : +proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs
## source     : C:/Users/Arthur/Documents/rscripts/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 (red, green, blue and near infrared).

tr_r = brick("../Data/s2/True color.tiff")
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded ellps WGS 84 in CRS definition: +proj=merc +a=6378137
## +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
## +wktext +no_defs +type=crs
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum World Geodetic System 1984 in CRS definition
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)
## crs        : +proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs
## source     : C:/Users/Arthur/Documents/rscripts/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")
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded ellps WGS 84 in CRS definition: +proj=merc +a=6378137
## +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
## +wktext +no_defs +type=crs

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum World Geodetic System 1984 in CRS definition
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")

It’s not yet possible to add composite rasters in leaflet (such as RGB). To include true color raster, we need to add it by using the mapview package and the viewRGB function. It does not allow to work the layer control in the same way.

viewRGB(tr_r, 1,2,3, map=m1, quantiles = c(0.05, 0.95))
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded ellps WGS 84 in CRS definition: +proj=merc +a=6378137
## +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
## +wktext +no_defs +type=crs
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum World Geodetic System 1984 in CRS definition
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded ellps WGS 84 in CRS definition: +proj=merc +a=6378137
## +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
## +wktext +no_defs +type=crs
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum World Geodetic System 1984 in CRS definition
## Warning in raster::projectRaster(x, raster::projectExtent(x, crs =
## sp::CRS(epsg3857)), : input and ouput crs are the same
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded ellps WGS 84 in CRS definition: +proj=merc +a=6378137
## +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
## +wktext +no_defs +type=crs
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum World Geodetic System 1984 in CRS definition

It is also possible to visualize multiple single band rasters at once using the leafsync package

m1 = leaflet(dfs) %>% 
  addTiles() %>% 
  addRasterImage(tr_r[[1]])
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded ellps WGS 84 in CRS definition: +proj=merc +a=6378137
## +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
## +wktext +no_defs +type=crs
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum World Geodetic System 1984 in CRS definition
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded ellps WGS 84 in CRS definition: +proj=merc +a=6378137
## +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
## +wktext +no_defs +type=crs
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum World Geodetic System 1984 in CRS definition
m2 = leaflet(dfs) %>%
  addTiles() %>% 
  addRasterImage(tr_r[[4]])
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded ellps WGS 84 in CRS definition: +proj=merc +a=6378137
## +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
## +wktext +no_defs +type=crs

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum World Geodetic System 1984 in CRS definition
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded ellps WGS 84 in CRS definition: +proj=merc +a=6378137
## +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
## +wktext +no_defs +type=crs
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum World Geodetic System 1984 in CRS definition
sync(m1,m2)

When multiple bands are available, it becomes possible to easily perform raster maths, allowing to obtain different indexes, such as the normalized difference vegetation index (NDVI). This index favors the visualisation of vegetation.

tr_r$NDVI=(tr_r$B08-tr_r$True_color.1)/(tr_r$B08+tr_r$True_color.1)
plot(tr_r$NDVI)

It is also possible to limit the extent of the satellite data to our zone of interest as delimited by our campus polygon.

crop_tr=mask(tr_r,campus)
plot(crop_tr$NDVI)

From there, we can visualize the campus NDVI values back into the leaflet environment.

val = as.numeric(c(0:1))
pal = colorNumeric(c("red","yellow","green"), val, na.color="transparent")
                 
ndvi_map=leaflet(df) %>%
  addTiles() %>%
  addRasterImage(crop_tr$NDVI, colors = pal , opacity = 0.3) %>%
  addLegend(pal = pal, values = val, title = "NDVI") %>%
  addCircles(color="black",stroke=T, opacity = 1)
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded ellps WGS 84 in CRS definition: +proj=merc +a=6378137
## +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
## +wktext +no_defs +type=crs
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum World Geodetic System 1984 in CRS definition
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded ellps WGS 84 in CRS definition: +proj=merc +a=6378137
## +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
## +wktext +no_defs +type=crs
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum World Geodetic System 1984 in CRS definition
## Warning in colors(.): Some values were outside the color scale and will be
## treated as NA
## Assuming "lon" and "lat" are longitude and latitude, respectively
ndvi_map

Since we have many data points scattered around campus that already have a lot of information, it can be interesting to add some remote sensing data to our data frame. That allows us to look at many types of interesting relations in the data.

Step 5. Extracting data from rasters

mean_ndvi = extract(crop_tr$NDVI, dfs, buffer=10, small=T, fun=mean)
sd_ndvi = extract(crop_tr$NDVI, dfs, buffer=10, small=T, fun=sd)

dfs$ndvi=mean_ndvi
dfs$ndvi_sd=sd_ndvi

plot(dfs$ndvi~dfs$aerial_drymass)

plot(dfs$ndvi~dfs$est_age)

boxplot(dfs$ndvi~dfs$dom_sp)

boxplot(dfs$ndvi~dfs$zone_eco)

Now that we know how to display and manipulate spatial data, let’s put it into a nice shareable format and save it.

To do so, let’s redefine our colors, add layer control and some additional information when clicking on points, then let’s save the map into an html object.

Step 6. Sharing your awesome map

ndvi_val = as.numeric(c(-1:1))
ndvi_pal = colorNumeric("RdYlGn",
                        ndvi_val,
                        na.color = "transparent")

circ_pal = colorNumeric("viridis",
                        dfs$total_drymass)

map = leaflet(dfs) %>% 
  addScaleBar("bottomleft") %>% 
  addTiles(group = "BaseOSM") %>%
  addProviderTiles("Stamen.Toner",
                   group = "BaseStamen") %>% 
  addRasterImage(crop_tr$NDVI,
                 colors = ndvi_pal,
                 opacity = 0.6,
                 group = "NDVI") %>%
  addLegend("bottomright",
            pal = ndvi_pal,
            values = ndvi_val,
            title = "NDVI",
            group = "NDVI") %>% 
  addCircles(color = "black",
             opacity = 0.5,
             fillColor = ~circ_pal(tot_drymass),
             fillOpacity = 0.8,
             radius = dfs$n_trees,
             group = "drymass",
             stroke = F,
             popup = paste("n Trees:", dfs$n_trees, "<br>",
                           "mean DHP:", dfs$mean_dhp)) %>% 
  addLegend("bottomright",
            pal = circ_pal,
            values = ~tot_drymass,
            title = "drymass (kg / 0.04 ha)",
            group ="drymass") %>% 
  addLayersControl("bottomleft",
    baseGroups = c("BaseOSM","BaseStamen"),
    overlayGroups = c("NDVI","drymass")
  )
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded ellps WGS 84 in CRS definition: +proj=merc +a=6378137
## +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
## +wktext +no_defs +type=crs
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum World Geodetic System 1984 in CRS definition
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded ellps WGS 84 in CRS definition: +proj=merc +a=6378137
## +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
## +wktext +no_defs +type=crs
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum World Geodetic System 1984 in CRS definition
map
saveWidget(map, file="../Figures/map.html", selfcontained = T)

Step 7. Sharing spatial data

There are many ways to share spatial data without loosing its spatial properties (CRS, for example).

The classics are the ESRI shapefiles, but they are very clunky, slow and require multiple files. Alternatives include the GeoJSON format, and more recently, geopackage (GPKG, which is also faster).

st_write(dfs, "../Data/out_df.gpkg","campus_trees",driver="GPKG", overwrite_layer = T, append=F)
## Deleting layer `campus_trees' using driver `GPKG'
## Writing layer `campus_trees' to data source `../Data/out_df.gpkg' using driver `GPKG'
## Writing 36 features with 18 fields and geometry type Point.
test = st_read("../Data/out_df.gpkg")
## Reading layer `campus_trees' from data source `C:\Users\Arthur\Documents\rscripts\NumeriLab_RGIS\Data\out_df.gpkg' using driver `GPKG'
## Simple feature collection with 36 features and 18 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -72.58601 ymin: 46.34263 xmax: -72.57278 ymax: 46.35089
## Geodetic CRS:  WGS 84
test
## Simple feature collection with 36 features and 18 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -72.58601 ymin: 46.34263 xmax: -72.57278 ymax: 46.35089
## Geodetic CRS:  WGS 84
## First 10 features:
##    zone parcelle             zone_eco sup_zone_ha sup_parc_ha tot_drymass
## 1    18      18A jeune pinede blanche       0.213        0.04        8861
## 2     4       4A jeune pinede blanche       0.552        0.04       11379
## 3    13      13A                mixte       1.441        0.04        6035
## 4    17      17A                mixte       2.789        0.04       13721
## 5    17      17B                mixte       2.782        0.04        9865
## 6    17      17C                mixte       2.782        0.04       10449
## 7    17      17D                mixte       2.782        0.04        6449
## 8    20      20A        mixte suranne       1.761        0.04       13441
## 9    20      20B        mixte suranne       1.761        0.04       20148
## 10   20      20C        mixte suranne       1.761        0.04       13799
##    aerial_drymass root_drymass n_trees total_C n_living_trees mean_dhp max_dhp
## 1            6877         1984      95    4735             84     14.9    29.0
## 2            8896         2483      50    6236             50     23.5    33.9
## 3            4692         1343      51    3146             38     16.6    28.9
## 4           10766         2955      57    7651             42     22.2    79.1
## 5            7726         2140      40    5397             29     23.5    42.4
## 6            8200         2250      39    5727             34     22.8    49.2
## 7            5023         1426      49    3425             41     17.2    35.2
## 8           10583         2858      28    7375             27     23.8    55.4
## 9           15939         4209      18   11114             17     39.9    72.2
## 10          10865         2934      33    7316             32     22.6    68.8
##    mean_height  dom_sp est_age      ndvi     ndvi_sd                       geom
## 1         14.1     PIB      45 0.6459874 0.026004052   POINT (-72.586 46.34558)
## 2         13.6     PIB      45 0.5877693 0.050616639 POINT (-72.57578 46.34992)
## 3         13.6     BOP      45 0.6476816 0.009793101 POINT (-72.58229 46.34743)
## 4         16.9 PIB/PIG      45 0.6389911 0.006086868   POINT (-72.583 46.34656)
## 5         18.3     PIB      45 0.6749924 0.011763637 POINT (-72.58397 46.34653)
## 6         15.1     PIB      45 0.6688315 0.010359851 POINT (-72.58511 46.34597)
## 7         13.0 PIB/ERR      45 0.6963485 0.015353180  POINT (-72.58475 46.3465)
## 8         17.1   Mixte      50 0.6944277 0.021046173 POINT (-72.58601 46.34511)
## 9         20.0   Mixte      70 0.6421419 0.014284640 POINT (-72.58511 46.34481)
## 10        16.5 Feuillu      55 0.6705704 0.011261256  POINT (-72.5843 46.34298)