Where are climatically similar locations?
This is not a straightforward question because a range of conditions contribute to climate including rainfall, temperature, and the seasonal and daily variation in these properties.
A key set of spatial data that contains useful information on these ecoclimatic properties are the BIOCLIM variables frequently used in species distributinon models.
We can download the data which contain the following variables.
library(tidyverse)
"data/bioclim/bioclim_var_description.csv" %>%
read_csv() %>%
knitr::kable()
var | desc |
---|---|
bio01 | Annual Mean Temperature |
bio02 | Mean Diurnal Range (Mean of monthly (max temp - min temp)) |
bio03 | Isothermality (BIO2/BIO7) (×100) |
bio04 | Temperature Seasonality (standard deviation ×100) |
bio05 | Max Temperature of Warmest Month |
bio06 | Min Temperature of Coldest Month |
bio07 | Temperature Annual Range (BIO5-BIO6) |
bio08 | Mean Temperature of Wettest Quarter |
bio09 | Mean Temperature of Driest Quarter |
bio10 | Mean Temperature of Warmest Quarter |
bio11 | Mean Temperature of Coldest Quarter |
bio12 | Annual Precipitation |
bio13 | Precipitation of Wettest Month |
bio14 | Precipitation of Driest Month |
bio15 | Precipitation Seasonality (Coefficient of Variation) |
bio16 | Precipitation of Wettest Quarter |
bio17 | Precipitation of Driest Quarter |
bio18 | Precipitation of Warmest Quarter |
bio19 | Precipitation of Coldest Quarter |
We can load this data into R and plot them for Australia.
library(sf)
library(tmap)
library(raster)
# rough Australia shapefile
data(World)
aus = World %>%
filter(name == "Australia") %>%
st_transform(crs=4326)
World = World %>% filter(continent != "Antarctica")
# load bioclim variables
bioclim = stack(list.files("data/bioclim/", pattern = ".tif", full.names=T))
bioclim = crop(bioclim, extent(-180,180,-50,50))
names(bioclim) = str_pad(str_extract(names(bioclim), pattern = "(?<=bio_)\\d+"), 2, "left", "0")
bioclim = bioclim[[order(gsub("X","", names(bioclim)))]]
names(bioclim) = gsub("X", "bio", names(bioclim))
# make plot theme
mytmtheme = tm_layout(legend.position = c("left", "bottom"), bg.color = "#37314a",
scale = 0.5,outer.bg.color = "#37314a", frame = F)
# bioclim
tm_shape(crop(bioclim, aus)) +
tm_raster(style = "cont", pal=rainbow(10)) +
tm_facets(free.scales = T) +
mytmtheme
Unfortunately, there are many variables to consider which makes visual comparisons of the similarity between locations difficult. Luckily, these variables are highly correlated with each other and we can extract the principle component and reduce the dimensionality and redundancy in the data.
We can see that the first princple component explains 83% of the variation in the data.
vals = values(bioclim)
vals = vals[complete.cases(vals), ]
pc = princomp(vals, scores = T)
tibble(Component =1:19, Variance = (pc$sdev^2)) %>%
mutate(Proportion = Variance/sum(Variance)) %>%
knitr::kable()
Component | Variance | Proportion |
---|---|---|
1 | 9.117532e+05 | 0.8344917 |
2 | 1.021896e+05 | 0.0935301 |
3 | 4.469786e+04 | 0.0409102 |
4 | 2.342313e+04 | 0.0214383 |
5 | 7.654384e+03 | 0.0070058 |
6 | 1.714997e+03 | 0.0015697 |
7 | 6.085929e+02 | 0.0005570 |
8 | 2.377121e+02 | 0.0002176 |
9 | 1.780822e+02 | 0.0001630 |
10 | 5.842379e+01 | 0.0000535 |
11 | 4.159151e+01 | 0.0000381 |
12 | 1.205113e+01 | 0.0000110 |
13 | 9.801420e+00 | 0.0000090 |
14 | 4.912088e+00 | 0.0000045 |
15 | 4.633543e-01 | 0.0000004 |
16 | 3.083914e-01 | 0.0000003 |
17 | 4.224400e-02 | 0.0000000 |
18 | 1.327750e-02 | 0.0000000 |
19 | 0.000000e+00 | 0.0000000 |
We can plot this axis of variation as 10 discrete quantiles for the world
pc1 = bioclim[[1]]
names(pc1) = "PC1"
pc1[!is.na(values(pc1))] = pc$scores[,1]
tm_shape(pc1) +
tm_raster(style="quantile", n = 8, stretch.palette = F, midpoint=0,
pal=RColorBrewer::brewer.pal(8, "Set2"))+
mytmtheme
The previous plot spans a huge range of conditions and we might want to focus on a particular location of interest such as Australia, Victoria, Melbourne. We can then compute the difference in PC1 scores between this location and other others. Here we truncate scores outside of those in the same decile of our location.
# melbourne loction and score
xy = matrix(c(144.946457, -37.840935), ncol=2)
z = raster::extract(pc1, xy)
diff = pc1 - z
legmax = quantile(abs(values(diff)), 0.1, na.rm=T)
diff[values(diff) > legmax] = NA
diff[values(diff) < -legmax] = NA
tm_shape(aus) +
tm_fill("grey") +
tm_shape(crop(diff, aus)) +
tm_raster(title = "PC1 (deviation from Melbourne)",
style="cont", n = 10, stretch.palette = F, midpoint=0,
pal=RColorBrewer::brewer.pal(10, "RdYlGn")) +
mytmtheme
And for the world.
tm_shape(World) +
tm_fill("grey") +
tm_shape(diff, is.master = T) +
tm_raster(title = "PC1 (deviation from Melbourne)",
style="cont", n = 10, stretch.palette = F, midpoint=0,
pal=RColorBrewer::brewer.pal(10, "RdYlGn")) +
mytmtheme