Climatic similarity

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

# bioclim
tm_shape(crop(bioclim, aus)) + 
  tm_raster(style = "cont", pal=rainbow(10)) + 
  tm_facets(free.scales = T) + 
  tm_layout(outer.bg.color = "black", )

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 components can explain 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 these plot this axis of variation as 10 discrete quantiles for the globe.

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"))+ 
  tm_layout(outer.bg.color = "black")

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, Grampians (Gariwerd). 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(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")) +
  tm_shape(aus) + 
  tm_borders() +
  tm_layout(outer.bg.color = "black")

And for the world.

tm_shape(diff) + 
  tm_raster(title = "PC1 (deviation from Melbourne)", 
            style="cont", n = 10, stretch.palette = F, midpoint=0,
            pal=RColorBrewer::brewer.pal(10, "RdYlGn")) + 
  tm_shape(World) + 
  tm_borders() +
  tm_layout(outer.bg.color = "black")

Related

comments powered by Disqus