遥感Python实战: 栅格数据逐像元偏相关分析(附代码)

导入上述源代码文件（核心偏相关分析文件）import pcorr pc = pcorr.ParCorr(data_dir=r'E:\公众号\偏相关分析\vars', min_years=5)

ParCorr_Res = pc.PCor(x='SOT',  y='GUD', covar=['ALT'],  method='pearson')

pc.plot_results(ParCorr_Res)

pc.save_results_to_tif(ParCorr_Res, r'E:\公众号\偏相关分析\结果\ParCorr_SOT_GUD')

# -*- coding: utf-8 -*-import osimport reimport numpy as npimport rasterioimport pandas as pdimport pingouin as pgimport warningsimport gcfrom rasterio.windows import Windowclass ParCorr(object):    """    Pixel-wise Partial Correlation (Optimized for Large Scale)    """    def __init__(self, data_dir, nodata=-99999, min_years=7):        self.data_dir = data_dir        self.nodata = nodata        self.min_years = min_years        self.profile = None  # 初始为 None    def _init_profile(self, var_name):        """        轻量级获取元数据：只打开第一个文件读取属性，不加载像素数据        """        var_path = os.path.join(self.data_dir, var_name)        if not os.path.exists(var_path):            raise FileNotFoundError(f"找不到文件夹: {var_path}")        tifs = sorted([f for f in os.listdir(var_path) if f.endswith(".tif")])        if not tifs:            raise FileNotFoundError(f"在 {var_path} 中没有找到 .tif 文件")        with rasterio.open(os.path.join(var_path, tifs[0])) as src:            self.profile = src.profile.copy()            self.rows = src.height            self.cols = src.width        print(f"✅ 元数据初始化成功: 范围 {self.rows}x{self.cols}")    def _load_var_window(self, var_path, window, pattern=r'\d{4}'):        """按窗口读取局部数据栈"""        files = sorted(            [f for f in os.listdir(var_path) if f.endswith(".tif")],            key=lambda x: int(re.findall(pattern, x)[0])        )        stack = []        for f in files:            with rasterio.open(os.path.join(var_path, f)) as src:                arr = src.read(1, window=window).astype(np.float32)                if src.nodata is not None:                    arr[arr == src.nodata] = np.nan                elif self.nodata is not None:                    arr[arr == self.nodata] = np.nan                stack.append(arr)        return np.stack(stack, axis=0)    def _pixel_partial_corr(self, x, y, covars, method):        """        x, y: (T,)        covars: dict{name: (T,)}        """        data = {'V1': x, 'V2': y}        for k, v in covars.items():            data[k] = v        df = pd.DataFrame(data)        df = df.dropna()        n = len(df)        if n < self.min_years:            return np.nan, np.nan, (np.nan, np.nan), np.nan        try:            res = pg.partial_corr(                data=df,                x='V1',                y='V2',                covar=list(covars.keys()),                method=method)            row = res.iloc[0]             r = float(row['r'])            p = float(row['p_val'])            ci_low, ci_high = row['CI95']            if not (-1 <= r <= 1):                return np.nan, np.nan, (np.nan, np.nan), np.nan            return n, r, (ci_low, ci_high), p        except Exception:            return np.nan, np.nan, (np.nan, np.nan), np.nan    def PCor(self, x, y, covar, method='pearson', pattern=r'\d{4}', chunk_size=50):        if isinstance(covar, str):            covar = [covar]        var_list = [x, y] + covar        if self.profile is None:            self._init_profile(x)        # 1. 预先获取元数据 (假设初始化时已通过某个文件获取了 self.profile)        # 如果没有，需要先打开一个文件获取行列数        rows = self.profile['height']        cols = self.profile['width']        # 2. 初始化结果地图 (使用 float32 节省 50% 内存)        res_maps = {            'n': np.full((rows, cols), np.nan, dtype=np.float32),            'r': np.full((rows, cols), np.nan, dtype=np.float32),            'p_val': np.full((rows, cols), np.nan, dtype=np.float32),            'CI95_low': np.full((rows, cols), np.nan, dtype=np.float32),            'CI95_high': np.full((rows, cols), np.nan, dtype=np.float32)        }        print(f"开始分块计算 (块大小: {chunk_size} 行): {rows}x{cols} 网格...")        # 3. 外层分块循环        for start_row in range(0, rows, chunk_size):            end_row = min(start_row + chunk_size, rows)            curr_height = end_row - start_row            window = Window(0, start_row, cols, curr_height)            stacks = {                v: self._load_var_window(os.path.join(self.data_dir, v), window, pattern)                for v in var_list            }            for r_offset in range(curr_height):                global_row = start_row + r_offset                for c0 in range(cols):                    x_ts = stacks[x][:, r_offset, c0]                    y_ts = stacks[y][:, r_offset, c0]                    cov_ts = {cv: stacks[cv][:, r_offset, c0] for cv in covar}                    # 基础空值检查 (保留原逻辑)                    if np.all(np.isnan(x_ts)) or np.all(np.isnan(y_ts)):                        continue                    if any(np.all(np.isnan(v)) for v in cov_ts.values()):                        continue                    try:                        with warnings.catch_warnings():                            warnings.simplefilter("error", RuntimeWarning)                            n, r, ci, p = self._pixel_partial_corr(x_ts, y_ts, cov_ts, method)                        if not np.isnan(r):                            res_maps['n'][global_row, c0] = n                            res_maps['r'][global_row, c0] = r                            res_maps['p_val'][global_row, c0] = p                            res_maps['CI95_low'][global_row, c0] = ci[0]                            res_maps['CI95_high'][global_row, c0] = ci[1]                    except RuntimeWarning as e:                        print("\n" + "="*50)                        print(f"❌ 数值错误 at Global Row: {global_row}, Col: {c0}")                        print(f"⚠️ 详情: {e}")                        print(f"🔍 X方差: {np.var(x_ts):.6f}, Y方差: {np.var(y_ts):.6f}")                        raise RuntimeError(f"计算终止") from e                    except Exception as e:                        raise e            del stacks            gc.collect()            print(f"进度: {end_row}/{rows} 行已完成", end='\r')        print("\n✅ 处理完成。")        return res_maps    def plot_results(self, results, figsize=(20, 12)):        """        绘制组图：将 n, r, p_val, CI95_low, CI95_high 组合展示        """        plot_config = {            'n': {'title': 'Sample Size (n)', 'cmap': 'viridis', 'vmin': self.min_years, 'vmax': None},            'r': {'title': 'Partial Corr (r)', 'cmap': 'RdBu_r', 'vmin': -1, 'vmax': 1},            'p_val': {'title': 'P-Value', 'cmap': 'YlGnBu_r', 'vmin': 0, 'vmax': 0.1},            'CI95_low': {'title': '95% CI Lower', 'cmap': 'RdBu_r', 'vmin': -1, 'vmax': 1},            'CI95_high': {'title': '95% CI Upper', 'cmap': 'RdBu_r', 'vmin': -1, 'vmax': 1}        }        keys = ['n', 'r', 'p_val', 'CI95_low', 'CI95_high']        fig, axes = plt.subplots(2, 3, figsize=figsize)        axes = axes.flatten()        for i, key in enumerate(keys):            data = results[key]            cfg = plot_config[key]            im = axes[i].imshow(data, cmap=cfg['cmap'],                                 vmin=cfg['vmin'], vmax=cfg['vmax'])            axes[i].set_title(cfg['title'], fontsize=14, fontweight='bold')            axes[i].axis('off') # 地理数据通常关闭坐标轴            plt.colorbar(im, ax=axes[i], orientation='vertical', fraction=0.046, pad=0.04)        if len(keys) < len(axes):            axes[-1].axis('off')        plt.tight_layout()        plt.show()    def save_results_to_tif(self, results, output_dir):        if not os.path.exists(output_dir):            os.makedirs(output_dir)        profile = self.profile.copy()        profile.update(            driver='GTiff',            dtype='float32',            count=1,            nodata=self.nodata,            compress='lzw')        for name, data in results.items():            out = data.astype(np.float32)            out[np.isnan(out)] = self.nodata            with rasterio.open(                os.path.join(output_dir, f'{name}.tif'),                'w',                **profile) as dst:                dst.write(out, 1)