Clarke van Steenderen

Clarke van Steenderen

Postdoctoral researcher in the Department of Zoology and Entomology, Centre for Biological Control, Rhodes University

Creating a climate map with GPS coordinate points

3 minute read

Published:

This R code shows you how to plot a map for a desired country or geographic region, and how to display a chosen climatic variable (e.g. mean annual temperature, annual precipitation) as an overlay. The WorldClim dataset has 19 variables to choose from.

The code below produces Figures A - D, where figures A and B show how to plot a map of South Africa, with average annual temperature and precipitation layered onto it. Some GPS coordinate points are added, just to show you how to add desired sampling sites if desired. Figures C and D show how you can filter the data, so that you only display regions that meet certain constraints. For example, perhaps you want to show only those areas where the mean annual temperature was less than 15 degrees Celsius, or areas where annual precipitation was less than 20 cm.

See the final maps at the end of this post 🗺️.


library(raster)
library(ggplot2)
library(janitor)
library(dplyr)

################################################
# Import GPS points from excel file on GitHub
################################################

sample_gps = readxl::read_xlsx("G:/My Drive/Tetramesa Project/Source Modifiers/source_modifiers_updated_10_06_2021.xlsx", sheet = "Tetramesa_Clarke")

# Apply some filtering specific to this datasheet
sample_gps = filter(sample_gps, mitochondrial_COI == "TRUE" | nuclear_28S == "TRUE")

# Check data import 
head(sample_gps)

# Clean the column names
sample_gps <- sample_gps %>%
  janitor::clean_names()

head(sample_gps)

################################################
# Access data from the WorldClim data set
################################################

climate = getData("worldclim", var="bio", res=2.5)

# Convert temperature metadata to degrees C, and precipitation from mm to cm
gain(climate)=0.1

# Plot a few maps to explore the data
plot(climate$bio1, main="Annual Mean Temperature")
plot(climate$bio5, main="Max Temperature of Warmest Month")
plot(climate$bio6, main="Min Temperature of Coldest Month")
plot(climate$bio3, main="Isothermality")

# Crop the world map to include just South Africa, and extract the mean annual temperature data ($bio1)
r1 <- crop(climate$bio1, extent(10,35,-35,-20))

# Extract the annual precipitation data ($bio12)
r2 = crop(climate$bio12, extent(10,35,-35,-20))

################################################################
# FIGURE A
################################################################

# Plot the mean annual temperature for South Africa
plot(r1, 
     col = rev( colorRampPalette(c("orange", "blue"))( 255 ) ),  
     main="Annual Mean Temperature (C): South Africa",
     xlab = "Longitude",
     ylab = "Latitude"
     )

# Add GPS coordinates
points(sample_gps$longitude, sample_gps$latitude, bg = "black", col = "white", pch = 21, cex = 2)

################################################################
# FIGURE B
################################################################

# Plot annual precipitation for South Africa
plot(r2, 
     col = rev( colorRampPalette(c("blue", "grey"))( 255 ) ),  
     main="Annual Precipitation (cm): South Africa",
     xlab = "Longitude",
     ylab = "Latitude"
)

# Add GPS coordinates
points(sample_gps$longitude, sample_gps$latitude, bg = "red", col = "black", pch = 21, cex = 2)

################################################################
# Experiment with filtering the data
################################################################

# Find the range of the mean annual temperatures
cellStats(r1,range)

# Find the range of annual precipitation
cellStats(r2,range)

# Select all the data where temperature is less than 15 degrees Celsius
range.temp = r1 < 15

# Select all the data where precipitation is less than 20 cm
range.ppt = r2 < 20

################################################################
# FIGURE C
################################################################

# Plot the filtered map: temp
plot(range.temp, legend = F, col = c("lightgrey", "blue"),
     main = "Annual mean temperature < 15 degrees C")

################################################################
# FIGURE D
################################################################

# Plot the filtered map: ppt
plot(range.ppt, legend = F, col = c("lightgrey", "red"),
     main = "Annual Precipitation < 20 cm")

# add some GPS points
points(sample_gps$longitude, sample_gps$latitude, bg = "black", col = "black", pch = 21, cex = 2)

Multiple sequence alignment Fasta file manipulation

6 minute read

Published:

This bit of code can be a great help to subset a Fasta file alignment based on certain conditions taken from a sample information file (e.g. subset only the samples collected in Spring, or those in Winter and greater than y cm, or everything except those in Autumn, etc.).

Plot a phylogeny in R using ggtree

2 minute read

Published:

This post will show you how to generate a phylogeny plot in R, with a heatmap alongside it. This example uses binary data, but it can be adapted to be plotted with a continuous variable.

Creating an elevation map with GPS coordinate points

1 minute read

Published:

This R script plots the map of a desired country, and colours it by elevation. It then adds GPS coordinate points to it (for example sampling sites) from an Excel file. Check out the final image at the bottom of this post 😃. The North arrow and Tetramesa wasp was added manually later in Inkscape.

Renaming all the sequences in a FASTA file automatically

2 minute read

Published:

I came across the problem of renaming sequences in a FASTA sequence alignment after downloading over 200 sequences from GenBank for four different genes. The sequence names were assigned by the GenBank accession number (e.g. MK1526) only. I wanted the sequence names to have the species name in it as well, for example MK1526_Canis_lupis. To avoid the tedious task of manually doing this, I wrote a few lines of R code that I hope will be of use to others with the same issue.