Global ecoclimatic similarity index

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()

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()

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