Skip to content

Latest commit

 

History

History
347 lines (267 loc) · 13.8 KB

mapping.md

File metadata and controls

347 lines (267 loc) · 13.8 KB
title author date output
Some mapping tools for environment layers
February 18-19, 2018, updated 2021-04-09
html_document
keep_md toc toc_depth
true
true
3

You are here: R workshop >> Session 6 Mapping >> mapping demo

GBIF data access using R

Session 6: Mapping demo

This session focuses on working with environment layers, mapping, cropping and masking layers using the Raster R-package and other tools.

GBIF data for taxon liverleaf (blåveis:no) - from Norway

library('rgbif') # rOpenSci r-package for GBIF data
library('mapr') # rOpenSci r-package for mapping (occurrence data)
sp_name <- "Hepatica nobilis"; kingdom <- "Plantae" # liverleaf (blaaveis:no), taxonKey=5371699
key <- rgbif::name_backbone(name=sp_name, kingdom=kingdom)$speciesKey
sp <- rgbif::occ_search(taxonKey=key, return="data", hasCoordinate=TRUE, country="NO", limit=100)
sp_m <- sp[c("name", "catalogNumber", "decimalLongitude","decimalLatitude", "basisOfRecord", "year", "municipality", "taxonKey", "occurrenceID")] ## Subset columns (for more useful map pop-up)
#rgbif::gbifmap(sp, region = "norway")
mapr::map_leaflet(sp_m, "decimalLongitude", "decimalLatitude", size=2, color="blue")

Map GBIF data for Hepatica nobilis from Norway

Read example species occurrence data into R

Or choose your own species data following the examples provided in session 3.

sp <- utils::read.delim("./demo_data/sp.txt", header=TRUE, dec=".", stringsAsFactors=FALSE) ## read from file
#head(sp, n=5) ## preview first 5 records
##utils::write.table(sp, file="./demo_data/sp.txt", sep="\t", row.names=FALSE, qmethod="double") ## write

Extract coordinates suitable for e.g. Maxent

library('utils')
xy <- sp[c("decimalLongitude","decimalLatitude")] ## Extract only the coordinates
sp_xy <- sp[c("species", "decimalLongitude","decimalLatitude")] ## Input format for Maxent
#utils::write.table(sp_xy, file="./demo_data/sp_xy.txt", sep="\t", row.names=FALSE, qmethod="double")
#readLines("./demo_data/sp_xy.txt", n=10)

Get (detailed) administrative borders for Norway from GADM

Slow to render, many details for the coast...

library(raster)
gadm_norway_1 <- getData('GADM', country='NOR', level=1, path="./demo_data") ## level 0,1,2,...
plot(gadm_norway_1, main="Adm. Boundaries Norway Level 1")
points(xy, col='blue', pch=20) ## plot species occurrence points to the map (smaller dots)
#legend("bottomright", title = "Legend", legend = "H. nobilis", pch = 20, pt.bg = "blue", bty = "n")

GADM Norway

Get (simpler) country borders for Norway from maptools

library(maptools)
library(rgeos)
data(wrld_simpl) ## vector/factor of ISO2, ISO2, NAME, REGION, SUBREGION, LON, LAT, ...
norway_mask <- subset(wrld_simpl, NAME=="Norway") ## extract the (simpler) border for Norway
plot(norway_mask, axes=FALSE, border="#777777") ## 
points(xy, col='blue', pch=20, cex=1) # plot species occurrence points to the map
title("Country mask for Norway from wrld_simpl")
legend("bottomright", title = "Legend", legend = "Occurrences", pch = 20, col="blue", cex = 0.9)

Border for Norway from maptool:wrld_simpl

Admin borders from DIVA-GIS (read shapefile)

DIVA-GIS data by country: [http://www.diva-gis.org/gdata]. DIVA-GIS admin borders for Norway: [http://biogeo.ucdavis.edu/data/diva/adm/NOR_adm.zip]. Notice that these are the very same (detailed) borders as available from GADM.

## Download DIVA-GIS admin borders for Norway
if (!file.exists('./demo_data/NOR_adm.zip')) {
    download.file('http://biogeo.ucdavis.edu/data/diva/adm/NOR_adm.zip', './demo_data/NOR_adm.zip')
    dir.create(file.path("./demo_data/NOR_adm")) ## create a folder for shape files
    unzip('./demo_data/NOR_adm.zip', exdir='./demo_data/NOR_adm')
}
library(raster)
library(sp)
NOR_adm0 <- raster::shapefile('./demo_data/NOR_adm/NOR_adm0.shp') ## country border
#NOR_adm1 <- raster::shapefile('./demo_data/NOR_adm/NOR_adm1.shp') ## county borders
#NOR_adm2 <- raster::shapefile('./demo_data/NOR_adm/NOR_adm2.shp') ## municipality borders
par(mfrow=c(1,3)) ## combining plots: nrows, ncols
##
plot(NOR_adm0, axes=FALSE, border="#666666")
points(xy, col='blue', pch=20, cex=1) ## pecies occurrence
title("NOR_adm0, country")
##
plot(NOR_adm1, axes=FALSE, border="#666666") ## may set line width -- lwd=INT, but increases rendering time substatially
points(xy, col='blue', pch=20, cex=1) ## species occurrence
title("NOR_adm1, county")
##
plot(NOR_adm2, axes=FALSE, border="#666666")
points(xy, col='blue', pch=20, cex=1) ## species occurrence
title("NOR_adm2, municipality")
#legend("bottomright", title = "Legend", legend = "Occurrences", pch = 20, col="blue", cex = 0.9)
##

DIVA-GIS NOR_adm (source: GADM)

Read environment layer from WorldClim into R

For more information about WorldClim see session 5.

require(raster) # spatial raster data management, works well with dismo
env <- getData('worldclim', var='bio', res=10) # 10 degree grid (approx 18.5 km, 342 km2 at equator) 85 MByte

Plot environment layers and species occurrences on a map

plot(env, 1, main=NULL, axes=FALSE) ## could add title here with main="Title"
title(main = bquote(italic(.(sp_name)) ~occurrences~on~Annual~mean~temperature~'(dCx10)'))
#plot(gadm_norway, add=TRUE) ## add admin county borders
points(xy, col='blue', pch=20) # plot species occurrence points to the map (smaller dots)

Bioclim 1, Annual mean temperature

plot(env, 12, main=NULL, axes=FALSE) # plot bioclim layer, BIO12 = Annual Precipitation
title(main = bquote(italic(.(sp_name)) ~occurrences~plotted~on~Annual~precipitation~'(mm)'))
axis(side=2, tck = -0.04, col.ticks="gray") ## add axis only left
points(xy, col='blue') # plot species occurrence points to the map

Bioclim 12, Annual precepitation

Crop and mask raster layers

Cutting large (global) environment layers to a smaller extent can save significant memory. If your species occurrence data are limited to a region (e.g. Norway, Scandinavia or similar) you might reduce computation time significantly by cropping your environment layers appropriatly.

Cut environment layer(s) to extent (result is always a square)

library(raster)
library(rgdal)
ext <- extent(3,35,54,72) ## minLon=3, maxLon=35, minLat=54, MaxLat=72 for Scandinavia
env_cut <- crop(env, ext, snap="out") ## square output
#env_cut <- crop(env, ext, snap="out", filename="./demo_data/env_cut.tif") ## add filename to save result
plot(env_cut)

Bioclim layers cropped to Scandinavia

Cut environment layer(s) to a mask (from a shapefile or other vector data)

data(wrld_simpl) ## here you can also read in your OWN vector data (e.g. study area)
norway_mask <- subset(wrld_simpl, NAME=="Norway")
env_crop <- crop(env, extent(norway_mask)) ## crop gives a square (cut to the extent of the mask)
env_mask <- mask(env_crop, norway_mask) ## mask removes data outside the mask
plot(env_mask)
plot(norway_mask, add=TRUE, lwd=2)

Bioclim layers masked to Norway

Plot with extent Scandinavia (using zoom)

Using zoom, the raster data in R workspace environment is still the same size. You only zoom into the region of interest for more useful maps.

library(raster)
ext <- extent(3,35,54,72) ## minLon=3, maxLon=35, minLat=54, MaxLat=72 for Scandinavia
plot(zoom(env$bio12, ext), add = TRUE); title(main="Bio12, Annual precipitation") 
points(xy, col='blue', pch=20) ## plot species occurrence points # Error msg: plot.new not called yet

Bioclim 12, Annual precepitation

Map of the cropped and masked raster layers

par(mfrow=c(1,2)) ## combining two plot with par(n_rows, n_columns)
##plot(env_cut$bio12); title(main="Bio12 cut to Scandinavia")
plot(env_crop$bio12); title(main="Bio12 cropped to Norway")
points(xy, col='blue', pch=20) ## add species occurrence points to cropped map
plot(env_mask$bio12); title(main="Bio12 masked to Norway")
points(xy, col='blue', pch=20) ## add species occurrence points to masked map

Bioclim 12 (Annual precepitation) respectively cropped and masked to Norway

Base map from Google

library(dismo)
g_no <- gmap("Norway")
trd <- geocode('Trondheim, Norway')
osl <- geocode('Oslo, Norway')
trd_merc <- Mercator(trd[, c('longitude', 'latitude')])
osl_merc <- Mercator(osl[, c('longitude', 'latitude')])
xy_merc <- Mercator(xy)
##
plot(g_no, interpolate=TRUE)
points(xy_merc, pch=20, col='blue') ## add species occurrence
points(trd_merc, pch='*', col='red', cex=3) ## add * for Trondheim
points(osl_merc, pch='*', col='red', cex=3) ## add * for Oslo

Base map from Google, with H. nobilis, Trondheim, Oslo

Testing remote sensing image for Trondheim downloaded from Landsat

I used the USGS LandsatLook Viewer and Sentinel2Look Viewer to download sattelite data for Trondheim.

## A rasterBrick is a multi-layer raster object (similar to a rasterStack)
rs_l <- brick('./demo_data/landsat_trondheim_web_mercartor_wgs84.tif')
rs_s <- brick('./demo_data/sentinel_trondheim_web_mercartor_wgs84.tif')
nlayers(rs_l); nlayers(rs_s) ## result = 3 layers
crs(rs_l); crs(rs_s) ## +proj=merc +a=6378137 +b=6378137 ... +units=m ...
ncell(rs_l); ncell(rs_s) ## rs_l = 100738 ## rs_s = 204530
dim(rs_l); dim(rs_s) ## rs_l = 241 418   3 ## rs_s = 362 565   3
res(rs_l);res(rs_s) ## rs_l = 30.09113 30.09113 ## rs_s = 20.03305 20.03305
#par(mfrow=c(2,1)) ## combining two plot with par(n_rows, n_columns) ## if landsat & sentinel
#plotRGB(rs_l, stretch="lin", axes=FALSE, main="Landsat True Color Composite") ## Landsat example
plotRGB(rs_s, stretch="lin", axes=FALSE, main="Sentinel True Color Composite") ## Sentinel example

Remote sensing data, sentinel, Trondheim

Size of environment layer can be LARGE if using the finer resolutions

## object.size(env) ## read the space allocated in memory for an environment variable
## format(object.size(env), units = "auto") ## Auto reports multiples of 1024
## format(object.size(env), units = "auto", standard = "SI") ## SI use multiples of 1000
cat("Size of env =", format(object.size(env), units = "auto"))
cat("\nSize of env_cut =", format(object.size(env_cut), units = "auto"))
cat("\nSize of env_mask =", format(object.size(env_mask), units = "auto"))
cat("\nSize of gadm_norway_1 =", format(object.size(gadm_norway_1), units = "auto"))
#rm(env) ## save memory - especially useful if using finer resolutions
  • Size of env = 235.5 Kb
  • Size of env_cut = 13.7 Kb
  • Size of env_mask = 940.8 Kb
  • Size of gadm_norway_1 = 12.8 Mb

The BioClim layers:

  • BIO1 = Annual Mean Temperature (degree C multiplied by 10, to avoid decimal point)
  • BIO2 = Mean Diurnal Range (Mean of monthly (max temp – min temp))
  • BIO3 = Isothermality (BIO2/BIO7) (* 100)
  • BIO4 = Temperature Seasonality (standard deviation *100)
  • BIO5 = Max Temperature of Warmest Month
  • BIO6 = Min Temperature of Coldest Month
  • BIO7 = Temperature Annual Range (BIO5-BIO6)
  • BIO8 = Mean Temperature of Wettest Quarter
  • BIO9 = Mean Temperature of Driest Quarter
  • BIO10 = Mean Temperature of Warmest Quarter
  • BIO11 = Mean Temperature of Coldest Quarter
  • BIO12 = Annual Precipitation (mm)
  • BIO13 = Precipitation of Wettest Month (mm)
  • BIO14 = Precipitation of Driest Month (mm)
  • BIO15 = Precipitation Seasonality (Coe cient of Variation)
  • BIO16 = Precipitation of Wettest Quarter
  • BIO17 = Precipitation of Driest Quarter
  • BIO18 = Precipitation of Warmest Quarter
  • BIO19 = Precipitation of Coldest Quarter

GBIF data for taxon liverleaf (blåveis:no) - from Trondheim

library('rgbif') # rOpenSci r-package for GBIF data
library('mapr') # rOpenSci r-package for mapping (occurrence data)
sp_name <- "Hepatica nobilis"; kingdom <- "Plantae" # liverleaf (blaaveis:no), taxonKey=5371699
key <- rgbif::name_backbone(name=sp_name, kingdom=kingdom)$speciesKey
bb <- c(10.2,63.3,10.6,63.5) # Trondheim
#bb <- c(5.25, 60.3, 5.4, 60.4) # Bergen
#bb <- c(18.7, 69.6, 19.2, 69.8) # Tromsoe
#bb <- c(10.6, 59.9, 10.9, 60.0) # Oslo
sp_bb <- rgbif::occ_search(taxonKey=key, return="data", hasCoordinate=TRUE, country="NO", geometry=bb, limit=100)
sp_bb_m <- sp_bb[c("name", "catalogNumber", "decimalLongitude","decimalLatitude", "basisOfRecord", "year", "municipality", "taxonKey", "occurrenceID")] ## Subset columns
mapr::map_leaflet(sp_bb_m, "decimalLongitude", "decimalLatitude", size=2, color="blue")

Map GBIF data with bounding box for Trondheim

Navigate to workshop home (Github) or the GitHub.io pages.

Session 6 HTML version on GitHub.io here >> mapping.html