In hhwagner1/LandGenCourse: Interface for course "Landscape Genetic Data Analysis with R"

1. Overview {-}

This bonus vignette uses the packages sf, terra, and tmap, yet it can only scratch the surface. For a thorough introduction to geocomputation with R, see this excellent Gitbook: https://geocompr.robinlovelace.net/index.html.

a) Goals {-}

This bonus material expands the Worked Example to show:

• How to import and export ESRI shapefiles.
• Compatibility with packages sp and raster
• How to create static and interactive maps with tmap.
• How to plot a categorical raster with a predefined color scheme.

Try modifying the code to import your own data!

b) Required R packages {-}

library(LandGenCourse)
library(sf)
library(GeNetIt)
library(terra)
library(tmap)
library(dplyr)
library(tibble)
library(here)

2. Import and export ESRI shapefiles {-}

ESRI shapefiles are a widely used data format for sharing geospatial vector data. With the package sf, they are easy to import and export. Here' we will export Sites.sf to a shapefile and then import it again.

a) Export 'sf' object to shapefile

The following code may produce a warning that column names were abbreviated. It writes the component files for the ESRI shapefile into the pre-existing folder output (the first line will create it if does not exist yet). Remember to remove the hashtags '#' to un-comment the code before running it.

The argument delete_dsn specifies whether any existing file with the same name should be deleted first (i.e., overwritten).

data(ralu.site)
#if(!dir.exists(here("output"))) dir.create(here("output"))
#dir.create(here("output/Sites"))
#st_write(ralu.site, here("output/Sites/Sites.shp"), delete_dsn = TRUE)

Navigate to the Sites folder in the output folder (Files tab). You should now see four component files of the shapefile. These four files are required parts of the shapefile, always keep (or share) them together! Some shapefiles will have additional, optional component files.

• Sites.dbf
• Sites.prj
• Sites.shp
• Sites.shx

b) Import shapefile to 'sf' object

Importing a shapefile is very easy with the function st_read from the sf package. While you only need to provide the path to the .shp file, the other files listed above must be in the same folder.

#Sites.sf_a <- st_read(here("output/Sites/Sites.shp"))
#Sites.sf_a 

3. Compatibility with sp and raster objects {-}

The packages sf and terra recently replaced the older packages sp and raster. You may still encounter code that requires objects from these packages. Here we show how to convert between sf and sp, and between terra and raster. The conversion is easy, though please note that sp objects are S4 objects.

a) Converting between sf and sp

It is easy to convert an sf object into a Spatial object of package sp, using the function as_Spatial of the sf package. For point data, the resulting class will be a SpatialPointsDataFrame.

data(ralu.site)
Sites.sp <- sf::as_Spatial(ralu.site)
Sites.sp

We can convert the Spatial object back to an sf object with the function st_as_sf of the sf package:

Sites.sf_b <- sf::st_as_sf(Sites.sp)
Sites.sf_b

b) Converting between terra and raster

To convert a SpatRaster (terra) with multiple layers into RasterStack of the package raster, we use the function stack of the raster package:

RasterMaps <- rast(system.file("extdata/covariates.tif", package="GeNetIt"))

RasterMaps.r <- raster::stack(RasterMaps)
RasterMaps.r 

To convert a single layer into a RasterLayer of the package raster, we would use the function raster:

nlcd.r <- raster::raster(RasterMaps$nlcd)
nlcd.r 

To convert from a RasterLayer (package raster) to a SpatRaster (package terra), we use the function rast of the terra package:

nlcd <- terra::rast(nlcd.r)
nlcd

We can use the same function to convert from a RasterStack (package raster) to a SpatRaster with multiple layers (package terra):

RasterMaps_b <- terra::rast(RasterMaps.r)
RasterMaps_b

4. Plotting spatial data with tmap {-}

a) Plot geometry {-}

The package sf makes a clear distinction between the geometry information (spatial coordinates: where in space) and attribute information (what's at these locations). Hence, when using the function plot with sf objects, we need to decide what we want to plot: geometry or attributes?

To plot the geometry, we use function st_geometry to extract the geometry information from the sf object:

data(ralu.site)
Sites.sf_c <- ralu.site
plot(st_geometry(Sites.sf_c))

b) Plot attributes in space {-}

If we don't extract the geometry, then R will assume that we want to plot attribute data. The default is to plot the first ten attributes.

Here we set the point character pch to a filled circle, which is symbol #16. With cex=2, we define the symbol size.

For an overview of 'pch' symbol numbers, and colors, check: http://vis.supstat.com/2013/04/plotting-symbols-and-color-palettes/

par(mar=c(2,2,2,2))
plot(Sites.sf_c, pch=16, cex=2)

This is pretty cool! Let's have a closer look at two of the variables.

• Basin (left): this is a factor, and each factor level is assigned a different color.
• Depth_m (right): this is a quantitative variable, and R automatically uses a color ramp (from blue to pink to orange) to indicate variation in the values.

Note: To learn about options for the plot function for sf objects, access the help file by typing ?plot and select 'Plot sf object'.

plot(Sites.sf_c[,c("Basin", "Depth_m")], pch=16)

c) Create a static bubble plot with 'tmap' {-}

The tmap package (for plotting thematic maps) is a great tool for plotting maps. It is based on the grammar of graphics concepts, which take a bit of getting used to. Most importantly, we need the following parts:

• tmap_mode("plot"): to plot a static map (default)
• tm_shape: this function defines the data to be used.
• tm_sf: this function defines what information should be plotted and how.

Note: for the second part, there are many other functions for various types of data.

If we use tm_sf without arguments (i.e., with the default settings), we get a plot of the geometry:

tmap_mode("plot")
tm_shape(Sites.sf_c) + tm_sf()

We can indicate an attribute to plot it. Also, there is a special function tm_bubbles for bubble plots. Here we define the bubble size by wetland depth and bubble color by basin. In addition, we specify that the legend should be placed outside of the plot, on the right.

tmap_mode("plot")
tm_shape(Sites.sf_c) + tm_bubbles(size="Depth_m", col="Basin") +
  tm_layout(legend.outside=TRUE, legend.outside.position="right") 

Let's make the boundary box (map extent) a little larger so that the symbols are not cut off. First we extract the boundary box of Sites.sf_c and save it as Bbox.

Bbox = st_bbox(Sites.sf_c)
Bbox

Then we define the range along x and y coordinates (delta.x, delta.y), set a zoom factor (Zoom) and add that fraction of the range on each side. Unfortunately, there is no dedicated function for this so we do this manually:

delta.x <- Bbox[3] - Bbox[1]
delta.y <- Bbox[4] - Bbox[2]
Zoom <- 0.1
Bbox2 <- Bbox + c(-delta.x, -delta.y, delta.x, delta.y) * Zoom
Bbox2

Now we add the boundary box information as an argument bbox in function tm_shape.

Note that we write the figure into an object, Map1, then plot the Map1. This will help e.g. with exporting the map, or we can later add more layers to this map object with +.

tmap_mode("plot")
Map1 <- tm_shape(Sites.sf_c, bbox=Bbox2) + 
  tm_bubbles(size="Depth_m", col="Basin") +
  tm_layout(legend.outside=TRUE, legend.outside.position="right") 
Map1

d) Create an interactive bubble plot with 'tmap' {-}

Creating an interactive map with a basemap from the internet is not difficult. Un-comment the code below by removing the hashtag (#) and run it. Go ahead and play with the interactive map!

#tmap_mode("view")
#tm_shape(Sites.sf_c) + tm_bubbles(size="Depth_m", col="Basin") 

By default, R will include an interactive menu to toggle between "Esri.WorldGrayCanvas", "OpenStreetMap", and "Esri.WorldTopoMap". We can add more base maps to this selection by providing a list of servers with the function tm_basemap. The first one listed will be shown by default.

#tm_shape(Sites.sf_c) + tm_bubbles(size="Depth_m", col="Basin") +
#tm_basemap(server = c("Esri.WorldTopoMap", "Esri.WorldGrayCanvas", 
#                        "OpenStreetMap", "OpenTopoMap", "Esri.WorldImagery"))

Let's make a few more changes:

• We change the symbolization to fixed-size, filled circles with a black border, and plot them on top of filled circles that indicate the basin. Note that although both attributes are in the same dataset Sites.sf_c, we need to include a data statement for each attribute we add to the map (i.e., map layer), using tm_shape(Sites.sf).
• We include the argument bbox=Bbox2 for at least one of the map layers. This is not needed for the interactive map but will look better when exporting it as a static map.

Un-comment the code below by removing the hashtag (#) and run it. Go ahead and toggle between the base maps in the map below! Different base maps are suitable for different situations (data type, symbol type, size of study area) and purposes.

#tmap_mode("view")

#Map2 <- tm_shape(Sites.sf_c, bbox=Bbox2) +  tm_sf("Basin", size=2, border.col="black") +
#  tm_shape(Sites.sf_c) + tm_sf(size=0.8, col="Depth_m", 
#                             palette = "Blues", border.col="black") +
#  tm_basemap(server = c("Esri.WorldTopoMap", "Esri.WorldGrayCanvas", 
#                        "OpenStreetMap", "OpenTopoMap", "Esri.WorldImagery"))
#Map2

e) Export maps {-}

Save the map to folder output:

• Remember to un-comment the code by removing the hashtags (#).
• Interactive map saved to html file. This includes all the features of the interactive map Map2.
• Static map saved to png file, with height = 7 inches. This drops the base map and retains only the symbolized data.

Go ahead and check out the two files!

#if(!dir.exists(here("output"))) dir.create(here("output"))
#tmap_save(Map1, here::here("output/StaticMap.png"), height=7)
#tmap_save(Map2, here::here("output/InteractiveMap.html"))

Navigate to the folder output (Files tab) and check out the saved maps! For the dynamic map (Map2), select "open in a web browser". You can share this file with others, who can open it in their browser and interact with the map without access to R or your data!

5. Plot a categorical map with predefined color scheme {-}

a) Define the raster attribute table {-}

Now to a more tricky topic. Recall that the last raster layer in the Worked Example, nlcd, contains categorical land cover data that are coded numerically. The terra package actually misinterpreted them as numeric data.

Let's extract the categorical raster layer into a new object 'NLCD'. We can use as.factor to tell R that this is a categorical raster layer.

Here, we save the (numerical) raster layer nlcd as categorical raster (factor) NLCD.

RasterMaps <- rast(system.file("extdata/covariates.tif", package="GeNetIt"))
NLCD <- terra::as.factor(RasterMaps$nlcd)

What values occur in the raster? These are codes for cover types.

Here is a description of the cover types: https://www.mrlc.gov/data/legends/national-land-cover-database-2019-nlcd2019-legend

levels(NLCD)[[1]]

We will import a table with predefined colors (using hex color code) from the file Colortable_LULC.csv that is included with LandGenCourse.

This list has more entries (e.g., 21-24) than we need, because not all US land cover classes occur in the study area.

Check in the table below that the colors and cover types are stored as 'character'. (Luckily, since R 4.0, this is the new default for function read.csv). If they were coded as factors that could lead to errors later on.

ColTab <- read.csv(system.file("extdata", "Colortable_LULC.csv", 
                            package = "LandGenCourse"), header=TRUE)
ColTab

We join this list with the list of factor levels of NLCD to create a raster attribute table (RAT). An RAT is a table that contains attributes for each distinct value in a raster.

• Use function left_join to extract the corresponding row from ColTab for each row in Levels(NLCD)
• Indicate which columns should be used to match the rows: by=c("ID"="value"). Note that ID is of type numeric whereas nlcd is character. As value is numeric, we match it with ID.
• Copy the values from variable attribute to variable nlcd
• Select the variables that we will need.

RAT <- dplyr::left_join(levels(NLCD)[[1]], ColTab, by=c("ID"="value")) %>%
  mutate(nlcd=attribute) %>% 
  dplyr::select(ID, nlcd, color)

We replace the list of levels by the RAT:

levels(NLCD) <- RAT
NLCD

b) Plot the map with terra::plot

Now we can plot the map with the predefined color palette, using the plot function for SpatRaster objects:

plot(NLCD, col=RAT$color)
points(Sites.sf_c, pch=21, col="black", bg="white", cex=1)

c) Plot the map with tmap

• We first define the dataset with tm_shape
• With tm_raster, we tell tmap that we want to plot raster values
• style="cat" tells tmap to interpret values as categories
• palette defines the colors
• labels defines the labels to be used for the categories
• title defines the legend caption for the categories
• With tm_layout, we tell tmap where to place the legend (outside to the right of map)
• With tm_grid(lines=FALSE) we tell tmap to show the coordinates along x and y axes but to suppress grid lines that would be drawn on top of the map (change it to TRUE to see the effect)

Map3 <- tm_shape(NLCD) + 
  tm_raster(style="cat", palette=RAT$color, labels=RAT$nlcd,
            title="Land cover") +
  tm_layout(legend.outside=TRUE, legend.outside.position="right") 
Map3

Let's beef it up a bit. We can do so by adding layers to Map3:

• To add an additional layer with the sampling sites to the map, we define the data for the layer with tm_shape.
• We define how the points should be symbolized with tm_symbols.
• We add a North arrow with tm_compass , using default settings.
• We add a scale bar with tm_scale_bar and set the background to a semitransparent (bg.alpha=0.5) light gray (bg.color="lightgray").

Map4 <- Map3 + tm_shape(Sites.sf_c) +
  tm_symbols(size=0.4, col="yellow", border.col="red") +
  tm_compass() + tm_scale_bar(bg.color="lightgray", bg.alpha=0.5)
Map4

This time, we'll export the map to a PDF file.

#if(!dir.exists(here("output"))) dir.create(here("output"))
#tmap_save(Map4, here::here("output/RasterMap.pdf"), height=6, width=8)

# Detach all packages except for some basic ones:
LandGenCourse::detachAllPackages()

hhwagner1/LandGenCourse documentation built on Feb. 17, 2024, 4:42 p.m.