两个栅格数据叠加展示：10*10的二元色彩图

复现文章【https://doi.org/10.1038/s41586-022-04788-w】的10*10二元色彩图  ：

参考代码来自：用R绘制炫酷的二元色彩+山体阴影浮雕地图。这篇的代码处理的是矢量数据，并且做的是3*3的二元色彩图。进行了改进，以实现对两个栅格数据进行10*10的二元色彩展示。

代码如下：

一、读取栅格数据并转换为数据框

以全球的人口数据pop和一个全球的熵数据entr为例，读取数据后做预处理并分别转换为data.frame格式：

#加载包
library(rstudioapi) #
library(tidyverse) # ggplot2, dplyr, tidyr, readr, purrr, tibble
library(magrittr) # pipes# 包提供了一种称为"管道"（pipe）的操作，使代码更简洁和可读性更高。
library(lintr) # code linting  是一个用于代码静态分析和代码质量检查的包。
library(sf) # spatial data handling
library(raster) # raster handling (needed for relief)
library(viridis) # viridis color scale是一个用于生成美观的颜色映射的包。
library(cowplot) # stack ggplots  是一个用于组合和排列ggplot2图形的包
library(rmarkdown)

#读取两个栅格数据
popd2015<-rast("...文件路径.../popd2015.tif")
Entropy<-rast("...文件路径.../Entropy.tif")
#读取矢量边界数据用来裁剪两个栅格（这里直接加载世界地图矢量）
world.map <- ne_countries(returnclass = "sf") |>filter(continent != "Antarctica")
globalCountry <- vect(world.map) 

#预处理
popd2015_m <- trim(mask(popd2015, globalCountry))   #做mask
Entropy_m <- trim(mask(Entropy, globalCountry))

popd2015_r <- resample(popd2015_m, Entropy_m, method="bilinear")   #重采样
Entropy_m <- trim(mask(Entropy_m, popd2015_r))     #掩膜
Entropy_e <- extend(Entropy_m, popd2015_r)       #统一范围

二、合并为一个带坐标信息的数据框

#转换为带有坐标信息xy的data.frame格式
popd2015_df <- as.data.frame(popd2015_r, xy=TRUE)  
Entropy_df <- as.data.frame(Entropy_e, xy=TRUE)

#对数据进行归一化处理，统一度量方便展示：
# 定义归一化函数
normalize <- function(x) {
  return((x - min(x, na.rm = TRUE)) / (max(x, na.rm = TRUE) - min(x, na.rm = TRUE)))
}
#分别进行归一化处理-----
head(popd2015_df)   #查看数据列的名字为popd2015_30
popd2015_df$popd2015_30_normalized <- normalize(popd2015_df$popd2015_30)
head(Entropy_df)   #查看数据列的名字为sum
Entropy_df$sum_normalized <- normalize(Entropy_df$sum)

#按照xy的对应关系进行left_join
popentr_df <- left_join(popd2015_df, Entropy_df, by = c("x", "y"))
popentr_df <- na.omit(popentr_df)  #去掉na值
head(popentr_df)    #查看列名
names(popentr_df)<-c("x","y","pop","pop_n","entr","entr_n")   #重命名列名，这里xy是经纬度，pop_n和entr_n是归一化后的数据

三、设置分位数和10*10颜色变量

原理是给每一个像元按照分位数分组10*10，每一组分配一种颜色变量，因此，需要逐个设置每个位置的颜色代码。如果需要改变颜色，则需要改10*10个颜色代码。（颜色代码是把参考图给颜色识别网站，逐个取色识别的，然后复制粘贴了100次，这里直接把颜色代码给出来啦。）

同理需要4*4，5*5的，修改对应的length.out和颜色tibble就可以了

#设置10个分位数
quantiles_pop <- popentr_df %>%
  pull(pop_n) %>%
  quantile(probs = seq(0, 1, length.out = 11))

quantiles_entr <- popentr_df %>%
  pull(entr_n) %>%
  quantile(probs = seq(0, 1, length.out = 11))

#注意如果你的分位数中有重复的值，后面代码会报错可能需要处理
#举例子：给重复的加一个很小的值
quantiles_entr #检查分位数，看分位数的大小，根据分位数的大小拟定adjustment 
adjustment <- 0.0001   #以0.0001为例，如果数据更小，需要换成更小的数值
for(i in 2:length(quantiles_entr)) {
  if(quantiles_entr[i] == quantiles_entr[i-1]) {
    quantiles_entr[i] <- quantiles_entr[i] + adjustment
  }
  if(i>3 && quantiles_entr[i] == quantiles_entr[i-2]) {
    quantiles_entr[i] <- quantiles_entr[i-1] + adjustment  #防止连续三个分位数重复
  }
}
quantiles_entr#调整好后，检查分位数是否是递增的，如果不是，考虑改换adjustment或者改进循环的代码


#设置10*10颜色变量
bivariate_color_scale <- tibble(
  "1 - 10" ="#00FFFF","2 - 10" = "#11E6FD", "3 - 10" ="#23CDFB","4 - 10" = "#35B4FA","5 - 10" ="#479BF8","6 - 10"= "#5883F6","7 - 10" ="#6A6AF5","8 - 10" ="#7C51F3","9 - 10" ="#8E38F1","10 - 10" ="#A020F0",
  "1 - 9" = "#1BFDFB", "2 - 9" = "#2AE4F7", "3 - 9" = "#3ACCF4", "4 - 9" = "#4AB4F0","5 - 9" = "#5A9CED","6 - 9" = "#6A83E9","7 - 9" = "#7A6BE6","8 - 9" = "#8A53E2","9 - 9" = "#9A3BDF","10 - 9" = "#AA23DC",
  "1 - 8" = "#36FCF7", "2 - 8" = "#44E4F1", "3 - 8" = "#52CCEC", "4 - 8" = "#60B5E7","5 - 8" = "#6E9DE2","6 - 8" = "#7C85DC","7 - 8" = "#8A6ED7","8 - 8" = "#9856D2","9 - 8" = "#A63ECD","10 - 8" = "#B527C8",
  "1 - 7" = "#51FBF3", "2 - 7" = "#5DE3EC", "3 - 7" = "#69CCE5", "4 - 7" = "#75B5DE","5 - 7" = "#819ED7","6 - 7" = "#8E86D0","7 - 7" = "#9A6FC9","8 - 7" = "#A658C2","9 - 7" = "#B241BB","10 - 7" = "#BF2AB5",
  "1 - 6" = "#6CFAEF", "2 - 6" = "#76E3E6", "3 - 6" = "#80CCDD", "4 - 6" = "#8BB6D5","5 - 6" = "#959FCC","6 - 6" = "#A088C3","7 - 6" = "#AA72BB","8 - 6" = "#B55BB2","9 - 6" = "#BF44A9","10 - 6" = "#CA2EA1",
  "1 - 5" = "#88F9EB", "2 - 5" = "#90E2E0", "3 - 5" = "#98CCD6", "4 - 5" = "#A1B6CB","5 - 5" = "#A9A0C1","6 - 5" = "#B289B7","7 - 5" = "#BA73AD","8 - 5" = "#C35DA2","9 - 5" = "#CB4798","10 - 5" = "#D4318E",
  "1 - 4" = "#A3F8E7", "2 - 4" = "#A9E2DA", "3 - 4" = "#B0CCCE", "4 - 4" = "#B7B7C2","5 - 4" = "#BDA1B6","6 - 4" = "#C48BAA","7 - 4" = "#CB769E","8 - 4" = "#D16092","9 - 4" = "#D84A86","10 - 4" = "#DF357A",
  "1 - 3" = "#BEF7E3", "2 - 3" = "#C2E1D5", "3 - 3" = "#C7CCC7", "4 - 3" = "#CCB7B9","5 - 3" = "#D1A2AB","6 - 3" = "#D58C9E","7 - 3" = "#DA7790","8 - 3" = "#DF6282","9 - 3" = "#E44D74","10 - 3" = "#E93867",
  "1 - 2" = "#D9F6DF", "2 - 2" = "#DBE1CF", "3 - 2" = "#DFCCBF", "4 - 2" = "#E2B8B0","5 - 2" = "#E5A3A0","6 - 2" = "#E88E91","7 - 2" = "#EB7981","8 - 2" = "#EE6572","9 - 2" = "#F15062","10 - 2" = "#F43C53",
  "1 - 1" = "#F5F5DC", "2 - 1" = "#F6E0CA", "3 - 1" = "#F7CCB9", "4 - 1" = "#F8B8A8","5 - 1" = "#F9A496","6 - 1" = "#FA9085","7 - 1" = "#FB7C74","8 - 1" = "#FC6862","9 - 1" = "#FD5451","10 - 1" = "#FF4040",
  
) %>%
  gather("group", "fill")

四、“逐像元”分配颜色

使用cut函数将pop和entr的值分配到相应的分位数桶中（需要保证分位数唯一）

popentr_df1 <- popentr_df %>%
  mutate(
    pop_quantiles = cut(pop_n, breaks = quantiles_pop, include.lowest = TRUE, labels = FALSE),
    entr_quantiles = cut(entr_n, breaks = quantiles_entr, include.lowest = TRUE, labels = FALSE)
  ) %>%
  # 将分位数转换为数字，并创建一个新的组合列
  mutate(
    group = paste(as.numeric(pop_quantiles), "-", as.numeric(entr_quantiles))
    #group = paste(as.numeric(entr_quantiles), "-", as.numeric(pop_quantiles))
  ) %>%
  # 假设您已经有了一个颜色表bivariate_color_scale，其中包含group和对应的颜色值
  left_join(bivariate_color_scale, by = "group") # 将颜色值与主数据框合并

五、绘制地图和图例

#我这里用到了全球矢量边界
library(rnaturalearth)
coast <- ne_coastline(scale = "small", returnclass = "sf") %>% vect()
crs <- '+proj=longlat +datum=WGS84'

#绘制地图
map<-ggplot(popentr_df1) +
  geom_tile(aes(x = x, y = y,fill = fill)) +
  scale_fill_identity() +
  geom_spatvector(data=coast,fill=NA)+coord_sf(crs = crs,xlim=c(-160,165),ylim=c(-56,90))+
  #geom_sf(data = canton_geo, fill = NA, color = "white", size = 0.5) +    # 为行政区划绘制边界
  labs(x = NULL, y = NULL,) +
  #theme_minimal()+
  theme_bw()+
  theme(
    legend.position="none"
    )
print(map)


# bivariate_color_scale是颜色表，并且包含group和对应的颜色值
# 分离group列为pop和entr两个独立的列，并转换为整数类型
bivariate_color_scale <- bivariate_color_scale %>%
  separate(group, into = c("pop", "entr"), sep = " - ", convert = TRUE) %>%
  mutate(
    pop = as.integer(pop),
    entr = as.integer(entr)
  )

#绘制图例
legend <- ggplot(data = bivariate_color_scale) +
  geom_tile(mapping = aes(x = pop, y = entr, fill = fill)) +
  scale_fill_identity() +
  labs(x = "Higher population density ⟶️",
       y = "Higher entropy ⟶️") +
  theme_void() +  # 使用theme_void()来移除多余的元素
  theme(
    axis.title.x = element_text(size = 16, angle = 0),  # 旋转X轴标题
    axis.title.y = element_text(size = 16, angle = 90),   # 旋转Y轴标题
    plot.margin = margin(t = 10, r = 10, b = 30, l = 10)  # 调整图表边距
  ) +
  coord_fixed()  # 保持坐标比例固定
print(legend)

到这里就完成了，但地图和图例是分别绘制的，额外拼在一起就是这样的效果啦：