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的数据进行绘制):




目录
相关文章
|
20天前
|
数据采集 缓存 定位技术
网络延迟对Python爬虫速度的影响分析
网络延迟对Python爬虫速度的影响分析
|
2月前
|
数据采集 JSON 数据处理
抓取和分析JSON数据:使用Python构建数据处理管道
在大数据时代,电商网站如亚马逊、京东等成为数据采集的重要来源。本文介绍如何使用Python结合代理IP、多线程等技术,高效、隐秘地抓取并处理电商网站的JSON数据。通过爬虫代理服务,模拟真实用户行为,提升抓取效率和稳定性。示例代码展示了如何抓取亚马逊商品信息并进行解析。
抓取和分析JSON数据:使用Python构建数据处理管道
|
22天前
|
数据采集 存储 JSON
Python爬虫开发中的分析与方案制定
Python爬虫开发中的分析与方案制定
|
29天前
|
数据可视化 开发者 Python
Python GUI开发:Tkinter与PyQt的实战应用与对比分析
【10月更文挑战第26天】本文介绍了Python中两种常用的GUI工具包——Tkinter和PyQt。Tkinter内置于Python标准库,适合初学者快速上手,提供基本的GUI组件和方法。PyQt基于Qt库,功能强大且灵活,适用于创建复杂的GUI应用程序。通过实战示例和对比分析,帮助开发者选择合适的工具包以满足项目需求。
86 7
|
2月前
|
数据可视化 算法 Python
基于OpenFOAM和Python的流场动态模态分解:从数据提取到POD-DMD分析
本文介绍了如何利用Python脚本结合动态模态分解(DMD)技术,分析从OpenFOAM模拟中提取的二维切片数据,以深入理解流体动力学现象。通过PyVista库处理VTK格式的模拟数据,进行POD和DMD分析,揭示流场中的主要能量结构及动态特征。此方法为研究复杂流动系统提供了有力工具。
81 2
基于OpenFOAM和Python的流场动态模态分解:从数据提取到POD-DMD分析
|
29天前
|
存储 数据处理 Python
Python科学计算:NumPy与SciPy的高效数据处理与分析
【10月更文挑战第27天】在科学计算和数据分析领域,Python凭借简洁的语法和强大的库支持广受欢迎。NumPy和SciPy作为Python科学计算的两大基石,提供了高效的数据处理和分析工具。NumPy的核心功能是N维数组对象(ndarray),支持高效的大型数据集操作;SciPy则在此基础上提供了线性代数、信号处理、优化和统计分析等多种科学计算工具。结合使用NumPy和SciPy,可以显著提升数据处理和分析的效率,使Python成为科学计算和数据分析的首选语言。
32 3
|
29天前
|
存储 机器学习/深度学习 算法
Python科学计算:NumPy与SciPy的高效数据处理与分析
【10月更文挑战第26天】NumPy和SciPy是Python科学计算领域的两大核心库。NumPy提供高效的多维数组对象和丰富的数学函数,而SciPy则在此基础上提供了更多高级的科学计算功能,如数值积分、优化和统计等。两者结合使Python在科学计算中具有极高的效率和广泛的应用。
46 2
|
2月前
|
数据采集 机器学习/深度学习 搜索推荐
Python自动化:关键词密度分析与搜索引擎优化
Python自动化:关键词密度分析与搜索引擎优化
|
2月前
|
数据可视化 算法 JavaScript
基于图论的时间序列数据平稳性与连通性分析:利用图形、数学和 Python 揭示时间序列数据中的隐藏模式
本文探讨了如何利用图论分析时间序列数据的平稳性和连通性。通过将时间序列数据转换为图结构,计算片段间的相似性,并构建连通图,可以揭示数据中的隐藏模式。文章介绍了平稳性的概念,提出了基于图的平稳性度量,并展示了图分区在可视化平稳性中的应用。此外,还模拟了不同平稳性和非平稳性程度的信号,分析了图度量的变化,为时间序列数据分析提供了新视角。
61 0
基于图论的时间序列数据平稳性与连通性分析:利用图形、数学和 Python 揭示时间序列数据中的隐藏模式
|
2月前
|
自然语言处理 算法 数据挖掘
探讨如何利用Python中的NLP工具,从被动收集到主动分析文本数据的过程
【10月更文挑战第11天】本文介绍了自然语言处理(NLP)在文本分析中的应用,从被动收集到主动分析的过程。通过Python代码示例,详细展示了文本预处理、特征提取、情感分析和主题建模等关键技术,帮助读者理解如何有效利用NLP工具进行文本数据分析。
52 2