1 什么是偏相关分析
- 相关性分析用来反映要素之间的相关程度，以相关系数表示NDVI/EVI/NPP/NEP等与气象（气温和降水）因素的相关性，在简单相关分析的基础上固定某一要素，计算另外两个要素间的相关性，得出偏相关系数。大于0时正相关，小于0时负相关。
- 结合显著性（P值）可以得到：
显著正相关：相关系数>0 且 p值<0.05
不显著正相关：相关系数>0 且 p值≥0.05
不显著负相关：相关系数<0 且 p值≥0.05
显著负相关：相关系数<0 且 p值<0.05 - 偏相关分析公式如下：
![]()
2 代码实现逐像元偏相关分析
""" 逐像元偏相关分析，支持并行加速与显著性输出： - 输出4张栅格图（偏相关 + p值） https://mp.weixin.qq.com/s/dgZdeDpspc_nWfee09Ib5w """importosimportnumpyasnpimportrasteriofromsklearn.linear_modelimportLinearRegressionfromscipy.statsimportpearsonrfromtqdmimporttqdmfrommultiprocessingimportPool,cpu_count# ─── 函数：读取 TIF 文件───────────────────────────────────────defread_stack(folder):files=sorted([os.path.join(folder,f)forfinos.listdir(folder)iff.endswith('.tif')])stack=[]forpathintqdm(files,desc=os.path.basename(folder)):withrasterio.open(path)assrc:arr=src.read(1).astype(np.float32)nodata=src.nodataifnodataisnotNone:arr[arr==nodata]=np.nan stack.append(arr)stack=np.stack(stack)# (T, H, W)returnstack,src.meta.copy()# ─── 函数：偏相关计算（单像元） ─────────────────────────────defpartial_corr_with_p(x,y,z):x,y,z=x.reshape(-1,1),y.reshape(-1,1),z.reshape(-1,1)valid=np.isfinite(x[:,0])&np.isfinite(y[:,0])&np.isfinite(z[:,0])ifvalid.sum()<3:returnnp.nan,np.nan x_res=LinearRegression().fit(z[valid],x[valid]).predict(z[valid])y_res=LinearRegression().fit(z[valid],y[valid]).predict(z[valid])r,p=pearsonr((x[valid]-x_res).ravel(),(y[valid]-y_res).ravel())returnfloat(r),float(p)# ─── 函数：处理单行像元（用于并行） ─────────────────────────defprocess_one_row(args):i,NPP,TEM,PRE=args _,_,W=NPP.shape r_row_TEM=np.full(W,np.nan,dtype=np.float32)p_row_TEM=np.full(W,np.nan,dtype=np.float32)r_row_PRE=np.full(W,np.nan,dtype=np.float32)p_row_PRE=np.full(W,np.nan,dtype=np.float32)forjinrange(W):x=NPP[:,i,j]y1=TEM[:,i,j]y2=PRE[:,i,j]# 🌐 跳过该像元如果任一变量完全为空if(np.isnan(x).all()ornp.isnan(y1).all()ornp.isnan(y2).all()):continue# 默认值为 nanr1,p1=partial_corr_with_p(x,y1,y2)r2,p2=partial_corr_with_p(x,y2,y1)r_row_TEM[j]=r1 p_row_TEM[j]=p1 r_row_PRE[j]=r2 p_row_PRE[j]=p2returni,r_row_TEM,p_row_TEM,r_row_PRE,p_row_PRE# ─── 并行执行偏相关分析 ─────────────────────────────────────defpixelwise_partial_corr_parallel(NPP,TEM,PRE,workers=None):T,H,W=NPP.shape r_TEM=np.full((H,W),np.nan,dtype=np.float32)p_TEM=np.full((H,W),np.nan,dtype=np.float32)r_PRE=np.full((H,W),np.nan,dtype=np.float32)p_PRE=np.full((H,W),np.nan,dtype=np.float32)args=[(i,NPP,TEM,PRE)foriinrange(H)]withPool(processes=workersorcpu_count())aspool:forresultintqdm(pool.imap_unordered(process_one_row,args),total=H,desc='🔁 Parallel Rows'):i,r_row_TEM,p_row_TEM,r_row_PRE,p_row_PRE=result r_TEM[i,:]=r_row_TEM p_TEM[i,:]=p_row_TEM r_PRE[i,:]=r_row_PRE p_PRE[i,:]=p_row_PREreturn(r_TEM,p_TEM),(r_PRE,p_PRE)# ─── 函数：保存 GeoTIFF ─────────────────────────────────────defsave_raster(data,path,meta):meta.update(count=1,dtype='float32',nodata=np.nan)withrasterio.open(path,'w',**meta)asdst:dst.write(data,1)# ─── 主执行流程 ──────────────────────────────────────────────if__name__=='__main__':NDVI_dir=r"E:\AAWORK\20260310\data-裁剪02"TEM_dir=r"E:\AAWORK\20260310\tem-裁剪02"PRE_dir=r"E:\AAWORK\20260310\pre-裁剪02"result_dir=r"E:\AAWORK\20260310\results"print("🔍 正在加载数据...")NDVI,meta=read_stack(NDVI_dir)TEM,_=read_stack(TEM_dir)PRE,_=read_stack(PRE_dir)print("🧠 开始并行偏相关计算...")(r_TEM,p_TEM),(r_PRE,p_PRE)=pixelwise_partial_corr_parallel(NDVI,TEM,PRE,workers=12)print("💾 正在保存结果...")save_raster(r_TEM,os.path.join(result_dir,'partial_corr_NDVI_TEM.tif'),meta)save_raster(p_TEM,os.path.join(result_dir,'pvalue_NDVI_TEM.tif'),meta)save_raster(r_PRE,os.path.join(result_dir,'partial_corr_NDVI_PRE.tif'),meta)save_raster(p_PRE,os.path.join(result_dir,'pvalue_NDVI_PRE.tif'),meta)print("✅ 分析完成！输出保存在:",result_dir)
3 结合显著性（P值）得到四分类代码
# coding:utf-8importnumpyasnpimportpymannkendallasmkimportosfromosgeoimportgdal,gdalnumericdefread_tif(filepath):dataset=gdal.Open(filepath)col=dataset.RasterXSize# 图像长度row=dataset.RasterYSize# 图像宽度geotrans=dataset.GetGeoTransform()# 读取仿射变换proj=dataset.GetProjection()# 读取投影data=dataset.ReadAsArray()# 转为numpy格式data=data.astype(np.float32)no_data_value=data[0][0]# 设定NoData值data[data==no_data_value]=np.nan# 将NoData转换为NaNreturncol,row,geotrans,proj,datadefsave_tif(data,reference_file,output_file):ds=gdal.Open(reference_file)shape=data.shape driver=gdal.GetDriverByName("GTiff")dataset=driver.Create(output_file,shape[1],shape[0],1,gdal.GDT_Int16)# 保存的数据类型dataset.SetGeoTransform(ds.GetGeoTransform())dataset.SetProjection(ds.GetProjection())dataset.GetRasterBand(1).WriteArray(data)dataset.FlushCache()defread_tif02(file):data=gdalnumeric.LoadFile(file)data=data.astype(np.float32)a=data[0][0]data[data==a]=np.nanreturndata""" 分类方法说明 分类基于以下标准（显著性水平α=0.05）： 显著正相关：相关系数>0 且 p值<0.05 不显著正相关：相关系数>0 且 p值≥0.05 不显著负相关：相关系数<0 且 p值≥0.05 显著负相关：相关系数<0 且 p值<0.05 """if__name__=='__main__':# 输入文件夹路径partial_file=r"E:\AAWORK\20260310\results\partial_corr_NDVI_TEM.tif"pvalue_file=r"E:\AAWORK\20260310\results\pvalue_NDVI_TEM.tif"out_file=r"E:\AAWORK\20260310\results"col,row,geotrans,proj,partial_data=read_tif(partial_file)pvalue_data=read_tif02(pvalue_file)classified=np.zeros_like(partial_data,dtype=np.int16)classified[(partial_data>0)&(pvalue_data<0.05)]=1classified[(partial_data>0)&(pvalue_data>=0.05)]=2classified[(partial_data<0)&(pvalue_data>=0.05)]=3classified[(partial_data<0)&(pvalue_data<0.05)]=4# 输出文件路径save_tif(classified,partial_file,os.path.join(out_file,'partial_pvalue_corr_NDVI_TEM.tif'))print("Processing completed!")
