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()
|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|
|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()
We can these plot this axis of variation as 10 discrete quantiles for the globe.
pc1 = bioclim[] 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")