Python:如何基于滑动窗口进行气候因子间的相关系数分析?(逐像元)

简介: Python:如何基于滑动窗口进行气候因子间的相关系数分析?(逐像元)

最近处理一些气候因子的统计分析,遇到一些问题,记录一下。

01 常规的相关系数简单说明


在研究滑动窗口前,我们先来研究一下常规的相关系数分析,为了简化问题,假定我们现在有四川省范围内2019年每月的降水量数据和NDVI影像,我们想研究不同区域降水量和NDVI的相关性,即研究降水量和NDVI的空间关联性。


基本思路就是:对于每一个像元,都有1~12月共计12个值的降水量和NDVI指数,那么基于这两列数据我们可以计算出该像元位置降水量和NDVI的相关系数(这里我比较推荐使用斯皮尔曼系数,因为皮尔逊系数要求我们的这两列数据满足正态分布,但实际上我们往往都是假定我们的数据满足正态分布,这样计算的相关系数其实是有偏差的),对于每一个像元我们都可以计算出一个相关系数,如此便可以得到四川省范围内的相关系数的空间分布图。


思路讲完了,我们就可以得到如下的影像。


02 滑动窗口下的相关系数分析

重点来了,我们想看看在一定范围内的相关系数,怎么办?


说明:我当前有4个excel文件,分别为降水量数据,NDVI数据,地表温度数据,土壤水分数据。对于每一个excel文件,前5列表示地形因子数据,往后13列为2003年~2015年的降水量/NDVI/地表温度/土壤水分数据,每一行表示一个像元的地形因子值、气候因子值。


大致如下:


思路:假定滑动窗口为3*3(想想ArcGIS中空间分析工具-区域统计-焦点统计,有一点点类似),从第一个像元位置开始,对于该像元的移动窗口(包括其本身和周围8个像元),我们对于窗口内的每一个像元,都进行如下处理:获取窗口像元的降水量值(2003年~2015年共计13个值)和NDVI值,计算该像元位置的相关系数。那么对于该窗口我们一共计算了9个相关系数,然后我们计算这9个相关系数的平均值作为中心像元的相关系数值。类似的,对于降水量和地表温度、降水量和土壤水分均是如此操作。


以下是使用python实现的代码:

# @炒茄子  2023-05-18
import numpy as np
import pandas as pd
from scipy.stats import spearmanr
# 注意: 此处使用皮尔逊系数进行相关系数的计算, 但是皮尔逊系数要求数据满足正态分布, 但是实际上很多数据并不满足正态分布, 因此在此处尝试使用斯皮尔曼系数进行相关系数的计算
# 注意: 斯皮尔曼系数要求数据必须是有序的, 即数据具有可比较的大小关系
"""
当前程序用计算不同空间窗口下两两变量之间的相关系数,探索二者之间在空间上是否存在一定的关联。
"""
# 读取年平均降水量和 NDVI 数据、地表温度数据、土壤水分数据
precip_data_year = pd.read_excel(r'E:\PRCP\table\kinds_aver\precipitation_sc_year.xlsx', sheet_name='降水年平均值数据')
ndvi_data_year = pd.read_excel(r'E:\PRCP\table\kinds_aver\ndvi_sc_year.xlsx', sheet_name='NDVI年平均值数据')
temperature_data_year = pd.read_excel(r'E:\PRCP\table\kinds_aver\temperature_sc_year.xlsx', sheet_name='地表温度年平均值数据')
soil_moisture_data_year = pd.read_excel(r'E:\PRCP\table\kinds_aver\soil_moisture_sc_year.xlsx', sheet_name='土壤水分年平均值数据')
# 读取降水量数据和NDVI等数据(不读取地形因子和经纬度列)
_precip_data_year = precip_data_year.iloc[:, 5:]
_ndvi_data_year = ndvi_data_year.iloc[:, 5:]
_temperature_data_year = temperature_data_year.iloc[:, 5:]
_soil_moisture_data_year = soil_moisture_data_year.iloc[:, 5:]
# 创建空的三维数组,用于维度转换后的结果(83, 113, number_of_years)
precip_data_year = np.zeros((83, 113, len(_precip_data_year.columns)), dtype=np.float32)
ndvi_data_year = np.zeros((83, 113, len(_ndvi_data_year.columns)), dtype=np.float32)
temperature_data_year = np.zeros((83, 113, len(_temperature_data_year.columns)), dtype=np.float32)
soil_moisture_data_year = np.zeros((83, 113, len(_soil_moisture_data_year.columns)), dtype=np.float32)
# 将二维数组转换为三维数组
for i in range(len(_precip_data_year.columns)):
    precip_data_year[:, :, i] = _precip_data_year.iloc[:, i].values.reshape(83, 113)
    ndvi_data_year[:, :, i] = _ndvi_data_year.iloc[:, i].values.reshape(83, 113)
    temperature_data_year[:, :, i] = _temperature_data_year.iloc[:, i].values.reshape(83, 113)
    soil_moisture_data_year[:, :, i] = _soil_moisture_data_year.iloc[:, i].values.reshape(83, 113)
# 窗口大小列表
window_sizes = [5, 7, 9, 11, 13, 15, 17, 19, 21]
# 初始化相关系数DataFrame
correlation_dict = {}
for window_size in window_sizes:
    # 初始化相关系数矩阵(np.nan初始化)
    correlation_matrix = np.zeros((precip_data_year.shape[0], precip_data_year.shape[1], 3), dtype=np.float32)
    correlation_matrix[:] = np.nan
    # 对于每个窗口,计算斯皮尔曼相关系数
    for i in range(precip_data_year.shape[0] - window_size):
        for j in range(precip_data_year.shape[1] - window_size):
            # 提取窗口中的数据
            window_precipitation = precip_data_year[i:i + window_size, j:j + window_size, :]
            window_ndvi = ndvi_data_year[i:i + window_size, j:j + window_size, :]
            window_temperature = temperature_data_year[i:i + window_size, j:j + window_size, :]
            window_soil_moisture = soil_moisture_data_year[i:i + window_size, j:j + window_size, :]
            # 在窗口中的每个像素上计算相关系数
            corr_per_pixel = np.zeros((window_size, window_size, 3), dtype=np.float32)
            for k in range(window_size):
                for l in range(window_size):
                    corr_per_pixel[k, l, 0], _ = spearmanr(window_precipitation[k, l, :], window_ndvi[k, l, :])
                    corr_per_pixel[k, l, 1], _ = spearmanr(window_precipitation[k, l, :], window_temperature[k, l, :])
                    corr_per_pixel[k, l, 2], _ = spearmanr(window_precipitation[k, l, :], window_soil_moisture[k, l, :])
            # 计算窗口内的平均相关系数
            correlation_matrix[i, j, 0] = np.mean(corr_per_pixel[:, :, 0])
            correlation_matrix[i, j, 1] = np.mean(corr_per_pixel[:, :, 1])
            correlation_matrix[i, j, 2] = np.mean(corr_per_pixel[:, :, 2])
    # 将相关系数矩阵添加到correlation_matrices中
    correlation_dict['precipitation_ndvi_' + str(window_size)] = pd.Series(
        correlation_matrix[:, :, 0][~np.isnan(correlation_matrix[:, :, 0])].flatten())
    correlation_dict['precipitation_temperature_' + str(window_size)] = pd.Series(
        correlation_matrix[:, :, 1][~np.isnan(correlation_matrix[:, :, 1])].flatten())
    correlation_dict['precipitation_soil_moisture_' + str(window_size)] = pd.Series(
        correlation_matrix[:, :, 2][~np.isnan(correlation_matrix[:, :, 2])].flatten())
# 将相关系数矩阵保存到Excel文件中
df = pd.concat(correlation_dict, axis=1)
df.to_excel(r'E:\PRCP\table\window_corr\correlation_matrices.xlsx', index=False, sheet_name='年平均相关系数')

注意,建议运行时将窗口window_sizes的值进行编辑修改,由于我没有使用并行处理,所以如此多的窗口大小循环会导致运行时间长。

最后,我们将输出的excel数据进行小提琴图的绘制(只使用了其中窗口大小为5的数据进行绘制):




目录
相关文章
|
27天前
|
缓存 Rust 算法
从混沌到秩序:Python的依赖管理工具分析
Python 的依赖管理工具一直没有标准化,主要原因包括历史发展的随意性、社区的分散性、多样化的使用场景、向后兼容性的挑战、缺乏统一治理以及生态系统的快速变化。依赖管理工具用于处理项目中的依赖关系,确保不同环境下的依赖项一致性,避免软件故障和兼容性问题。常用的 Python 依赖管理工具如 pip、venv、pip-tools、Pipenv、Poetry 等各有优缺点,选择时需根据项目需求权衡。新工具如 uv 和 Pixi 在性能和功能上有所改进,值得考虑。
84 35
|
1月前
|
机器学习/深度学习 数据可视化 数据挖掘
使用Python实现基于矩阵分解的长期事件(MFLEs)时间序列分析
在现代数据分析中,高维时间序列数据的处理和预测极具挑战性。基于矩阵分解的长期事件(MFLEs)分析技术应运而生,通过降维和时间序列特性结合,有效应对大规模数据。MFLE利用矩阵分解提取潜在特征,降低计算复杂度,过滤噪声,并发现主要模式。相比传统方法如ARIMA和深度学习模型如LSTM,MFLE在多变量处理、计算效率和可解释性上更具优势。通过合理应用MFLE,可在物联网、金融等领域获得良好分析效果。
64 0
使用Python实现基于矩阵分解的长期事件(MFLEs)时间序列分析
|
28天前
|
数据采集 数据可视化 数据挖掘
金融波动率的多模型建模研究:GARCH族与HAR模型的Python实现与对比分析
本文探讨了金融资产波动率建模中的三种主流方法:GARCH、GJR-GARCH和HAR模型,基于SPY的实际交易数据进行实证分析。GARCH模型捕捉波动率聚类特征,GJR-GARCH引入杠杆效应,HAR整合多时间尺度波动率信息。通过Python实现模型估计与性能比较,展示了各模型在风险管理、衍生品定价等领域的应用优势。
251 66
金融波动率的多模型建模研究:GARCH族与HAR模型的Python实现与对比分析
|
18天前
|
并行计算 安全 Java
Python GIL(全局解释器锁)机制对多线程性能影响的深度分析
在Python开发中,GIL(全局解释器锁)一直备受关注。本文基于CPython解释器,探讨GIL的技术本质及其对程序性能的影响。GIL确保同一时刻只有一个线程执行代码,以保护内存管理的安全性,但也限制了多线程并行计算的效率。文章分析了GIL的必要性、局限性,并介绍了多进程、异步编程等替代方案。尽管Python 3.13计划移除GIL,但该特性至少要到2028年才会默认禁用,因此理解GIL仍至关重要。
97 16
Python GIL(全局解释器锁)机制对多线程性能影响的深度分析
|
1月前
|
数据可视化 算法 数据挖掘
Python时间序列分析工具Aeon使用指南
**Aeon** 是一个遵循 scikit-learn API 风格的开源 Python 库,专注于时间序列处理。它提供了分类、回归、聚类、预测建模和数据预处理等功能模块,支持多种算法和自定义距离度量。Aeon 活跃开发并持续更新至2024年,与 pandas 1.4.0 版本兼容,内置可视化工具,适合数据探索和基础分析任务。尽管在高级功能和性能优化方面有提升空间,但其简洁的 API 和完整的基础功能使其成为时间序列分析的有效工具。
80 37
Python时间序列分析工具Aeon使用指南
|
1月前
|
机器学习/深度学习 运维 数据可视化
Python时间序列分析:使用TSFresh进行自动化特征提取
TSFresh 是一个专门用于时间序列数据特征自动提取的框架,支持分类、回归和异常检测等机器学习任务。它通过自动化特征工程流程,处理数百个统计特征(如均值、方差、自相关性等),并通过假设检验筛选显著特征,提升分析效率。TSFresh 支持单变量和多变量时间序列数据,能够与 scikit-learn 等库无缝集成,适用于大规模时间序列数据的特征提取与模型训练。其工作流程包括数据格式转换、特征提取和选择,并提供可视化工具帮助理解特征分布及与目标变量的关系。
75 16
Python时间序列分析:使用TSFresh进行自动化特征提取
|
1月前
|
数据采集 缓存 API
python爬取Boss直聘,分析北京招聘市场
本文介绍了如何使用Python爬虫技术从Boss直聘平台上获取深圳地区的招聘数据,并进行数据分析,以帮助求职者更好地了解市场动态和职位需求。
|
2月前
|
机器学习/深度学习 数据采集 数据挖掘
使用Python实现智能食品消费市场分析的深度学习模型
使用Python实现智能食品消费市场分析的深度学习模型
153 36
|
2月前
|
数据可视化 算法 数据挖掘
Python量化投资实践:基于蒙特卡洛模拟的投资组合风险建模与分析
蒙特卡洛模拟是一种利用重复随机抽样解决确定性问题的计算方法,广泛应用于金融领域的不确定性建模和风险评估。本文介绍如何使用Python和EODHD API获取历史交易数据,通过模拟生成未来价格路径,分析投资风险与收益,包括VaR和CVaR计算,以辅助投资者制定合理决策。
117 15
|
2月前
|
机器学习/深度学习 数据采集 数据挖掘
使用Python实现智能食品消费趋势分析的深度学习模型
使用Python实现智能食品消费趋势分析的深度学习模型
156 18

热门文章

最新文章

推荐镜像

更多