Python遥感开发之偏相关分析

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):filessorted([os.path.join(folder,f)forfinos.listdir(folder)iff.endswith(.tif)])stack[]forpathintqdm(files,descos.path.basename(folder)):withrasterio.open(path)assrc:arrsrc.read(1).astype(np.float32)nodatasrc.nodataifnodataisnotNone:arr[arrnodata]np.nan stack.append(arr)stacknp.stack(stack)# (T, H, W)returnstack,src.meta.copy()# ─── 函数偏相关计算单像元 ─────────────────────────────defpartial_corr_with_p(x,y,z):x,y,zx.reshape(-1,1),y.reshape(-1,1),z.reshape(-1,1)validnp.isfinite(x[:,0])np.isfinite(y[:,0])np.isfinite(z[:,0])ifvalid.sum()3:returnnp.nan,np.nan x_resLinearRegression().fit(z[valid],x[valid]).predict(z[valid])y_resLinearRegression().fit(z[valid],y[valid]).predict(z[valid])r,ppearsonr((x[valid]-x_res).ravel(),(y[valid]-y_res).ravel())returnfloat(r),float(p)# ─── 函数处理单行像元用于并行 ─────────────────────────defprocess_one_row(args):i,NPP,TEM,PREargs _,_,WNPP.shape r_row_TEMnp.full(W,np.nan,dtypenp.float32)p_row_TEMnp.full(W,np.nan,dtypenp.float32)r_row_PREnp.full(W,np.nan,dtypenp.float32)p_row_PREnp.full(W,np.nan,dtypenp.float32)forjinrange(W):xNPP[:,i,j]y1TEM[:,i,j]y2PRE[:,i,j]# 跳过该像元如果任一变量完全为空if(np.isnan(x).all()ornp.isnan(y1).all()ornp.isnan(y2).all()):continue# 默认值为 nanr1,p1partial_corr_with_p(x,y1,y2)r2,p2partial_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,workersNone):T,H,WNPP.shape r_TEMnp.full((H,W),np.nan,dtypenp.float32)p_TEMnp.full((H,W),np.nan,dtypenp.float32)r_PREnp.full((H,W),np.nan,dtypenp.float32)p_PREnp.full((H,W),np.nan,dtypenp.float32)args[(i,NPP,TEM,PRE)foriinrange(H)]withPool(processesworkersorcpu_count())aspool:forresultintqdm(pool.imap_unordered(process_one_row,args),totalH,desc Parallel Rows):i,r_row_TEM,p_row_TEM,r_row_PRE,p_row_PREresult 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(count1,dtypefloat32,nodatanp.nan)withrasterio.open(path,w,**meta)asdst:dst.write(data,1)# ─── 主执行流程 ──────────────────────────────────────────────if__name____main__:NDVI_dirrE:\AAWORK\20260310\data-裁剪02TEM_dirrE:\AAWORK\20260310\tem-裁剪02PRE_dirrE:\AAWORK\20260310\pre-裁剪02result_dirrE:\AAWORK\20260310\resultsprint( 正在加载数据...)NDVI,metaread_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,workers12)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):datasetgdal.Open(filepath)coldataset.RasterXSize# 图像长度rowdataset.RasterYSize# 图像宽度geotransdataset.GetGeoTransform()# 读取仿射变换projdataset.GetProjection()# 读取投影datadataset.ReadAsArray()# 转为numpy格式datadata.astype(np.float32)no_data_valuedata[0][0]# 设定NoData值data[datano_data_value]np.nan# 将NoData转换为NaNreturncol,row,geotrans,proj,datadefsave_tif(data,reference_file,output_file):dsgdal.Open(reference_file)shapedata.shape drivergdal.GetDriverByName(GTiff)datasetdriver.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):datagdalnumeric.LoadFile(file)datadata.astype(np.float32)adata[0][0]data[dataa]np.nanreturndata 分类方法说明 分类基于以下标准显著性水平α0.05 显著正相关相关系数0 且 p值0.05 不显著正相关相关系数0 且 p值≥0.05 不显著负相关相关系数0 且 p值≥0.05 显著负相关相关系数0 且 p值0.05 if__name____main__:# 输入文件夹路径partial_filerE:\AAWORK\20260310\results\partial_corr_NDVI_TEM.tifpvalue_filerE:\AAWORK\20260310\results\pvalue_NDVI_TEM.tifout_filerE:\AAWORK\20260310\resultscol,row,geotrans,proj,partial_dataread_tif(partial_file)pvalue_dataread_tif02(pvalue_file)classifiednp.zeros_like(partial_data,dtypenp.int16)classified[(partial_data0)(pvalue_data0.05)]1classified[(partial_data0)(pvalue_data0.05)]2classified[(partial_data0)(pvalue_data0.05)]3classified[(partial_data0)(pvalue_data0.05)]4# 输出文件路径save_tif(classified,partial_file,os.path.join(out_file,partial_pvalue_corr_NDVI_TEM.tif))print(Processing completed!)