（6）计算两个栅格变量数据的Copula联合概率分布

---

1. 概述

前面我讲了大量关于二维和高维的Copula。有同学问我是否能在多维时间和空间上生成基于栅格尺度的联合概率分布。比如一个简单的问题，2000-2020年某个区域每日的高温和降水（干热）的联合概率分布，计算出整个区域的概率分布特征。

---

2. 基本原理和步骤

首先要学会读取矩阵数据，这里取决于你输入的数据类型是什么样的，可以是txt，csv，tiff等等，我们这里主要说明读取tif和nc格式的文件，大家可以去我的专栏看看如何读取处理tif格式、nc和hdf格式的文件，这是常见的带有空间信息的栅格文件。这种文件带有经纬度坐标、每个格子里有自己的数值，代表空间信息，再加上一维时间信息，就成了三维矩阵影像，看起来像一本厚厚的书，总页数（影像数量）取决于时间跨度和时间步长。
根据研究需求确定阈值和公式计算联合分布概率，公式我置顶的文章是有的。
这里以这种降水（X）低于某个阈值，温度（Y）高于某个阈值为例，生成联合概率分布，公式如下：

最后将结果输出成结果栅格矩阵，便于显示出图。

3. 代码出图

3.1 主要代码片段

% 清除所有变量和控制台输出，初始化执行环境
clear all
clc
tic; % 开始计时，用于测量脚本执行时间

% 读取GeoTIFF格式的历史平均降水量数据和最高温度数据
meanprehistoricdata = imread('H:\CMIP6\PRE\Historical\meanprehistoricdata.tif');
meantmaxhistoricdata = imread('H:\CMIP6\PRE\Historical\meantmaxhistoricdata.tif');

% 读取NetCDF格式的CanESM5模型历史降水量数据
pre = ncread('I:\CMIP6\PRE\Historical\CanESM5\BC_PREC_CanESM5_historical_1979.nc', 'pre');
tamx = ncread('I:\CMIP6\PRE\Historical\CanESM5\BC_PREC_CanESM5_historical_1979.nc', 'tamx');

% 获取meanprehistoricdata矩阵的尺寸
[rows, cols, depth] = size(meanprehistoricdata);

% 初始化一个NaN矩阵，用于存储干旱和高温的联合概率
dryhot = nan(rows, cols);

% 遍历每个像素点
for i = 1:rows
   for j = 1:cols
       predata = meanprehistoricdata(i, j, :); % 获取第i行第j列的所有降水量数据
       tmaxdata = meantmaxhistoricdata(i, j, :); % 获取第i行第j列的所有最高温度数据

       % 将数据转换为列向量
       pre = predata(:); % 降水量
       tmax = tmaxdata(:); % 最高温度

       % 去除异常值，只有当tmax和pre都不全是NaN时才处理
       if all(isnan(tmax)) || all(isnan(pre))
           P = nan;
       else
           % 筛选有效的数据点
           index = pre >= 0 & ~isnan(pre) & ~isnan(tmax);
           predata = pre(index);
           tmaxdata = tmax(index);

           % 使用边缘分布拟合数据
           [D_U1, PD_U1] = marginfitdist(predata);
           [D_U2, PD_U2] = marginfitdist(tmaxdata);

           % 计算拟合的边缘概率
           EP1 = cdf(PD_U1{1}, predata);
           spi = cdf(PD_U1{1}, -0.5); % 降水量阈值的累积分布函数值
           EP2 = cdf(PD_U2{1}, tmaxdata);
           sti = cdf(PD_U2{1}, 0.5); % 温度阈值的累积分布函数值

           EPpretmax = [EP1 EP2]; 
           EPpretmax(EPpretmax == 0) = 1e-4;
           [Family1, thetahat1, loglik1] = doublecopulaselect(EPpretmax, 'aic'); % 选择最佳的双变量copula模型

           % 计算干旱和高温的联合概率
           EP = [spi sti];
           P = spi - copulascdf(EP, Family1, thetahat1);
       end
       dryhot(i, j) = P; % 将计算结果存储
   end
end

% 如果需要，翻转矩阵以匹配地理坐标系
dryhot = flip(dryhot', 2);

% 读取经纬度数据
lon = ncread('I:\CMIP6\PRE\Historical\CanESM5\BC_PREC_CanESM5_historical_1979.nc', 'lon');
lat = ncread('I:\CMIP6\PRE\Historical\CanESM5\BC_PREC_CanESM5_historical_1979.nc', 'lat');

% 定义输出文件夹
outputfiles = 'I:\CMIP6\Results\';

% 定义输出的地理参考对象
R = georasterref('RasterSize', size(dryhot), 'Latlim', [double(min(lat)) double(max(lat))], 'Lonlim', [double(min(lon)) double(max(lon))]);

% 将结果写入GeoTIFF文件
geotiffwrite(fullfile(outputfiles, ['dryhot.tif']), dryhot, R);

disp('well done') % 输出执行完成的消息

toc; % 结束计时，输出脚本总执行时间

3.2 出图效果

欢迎关注gongzhonghao 趣品科研，回复 copula，获取更多代码和前沿论文资讯等相关内容

4. 参考文献

He, K., Chen, X., Zhou, J., Zhao, D., & Yu, X. (2024). Compound successive dry-hot and wet extremes in China with global warming and urbanization. Journal of Hydrology, 636, 131332
He, K., Chen, X., Xuan, Y., Chunyu, D., & Dongmei, Z. (2024). Evaluation and prediction of compound geohazards in highly urbanized regions across China’s Greater Bay Area. Journal of Cleaner Production, 141641, https://doi.org/10.1016/j.jclepro.2024.141641.
Liu, Z., T. Tornros, and L. Menzel (2016), A probabilistic prediction network for hydrological drought identification and environmental flow assessment, Water Resour. Res., 52, 6243–6262, doi:10.1002/2016WR019106.
吴海江. 基于气象干旱和高温的中国农业干旱预测模型研究[D].西北农林科技大学,2022.DOI:10.27409/d.cnki.gxbnu.2021.000534.
Jiang, T., Su, X., Zhang, G., Zhang, T., & Wu, H. (2023). Estimating propagation probability from meteorological to ecological droughts using a hybrid machine learning copula method. Hydrology and Earth System Sciences, 27(2), 559-576. https://doi.org/10.5194/hess-27-559-2023
赵晓阳. 黄河流域干热变化及对植被影响研究[D].河南大学,2023.DOI:10.27114/d.cnki.ghnau.2023.000182.
Hao, Z. et al., 2020. Impact of dependence changes on the likelihood of hot extremes under drought conditions in the United States. Journal of Hydrology, 581: 124410. DOI:https://doi.org/10.1016/j.jhydrol.2019.124410