Python遥感实战:基于多个自变量的岭回归分析

# -*- coding: utf-8 -*-# Packagesimport osimport reimport numpy as npimport rasterioimport pandas as pdfrom scipy import statsclass RidgeTrend(object):    """    Ridge Regression for Remote Sensing with Absolute and Relative Contribution analysis.    This class performs pixel-wise regression and exports results to GeoTIFF format.    """    def __init__(self, data_dir, nodata=-99999):        """        Initialize the RidgeTrend tool.        :param data_dir: Directory containing subfolders for each variable.        :param nodata: Value used for masking and output NoData.        """        self.data_dir = data_dir        self.nodata = nodata        self.profile = None    def _load_var_stack(self, var_path, pattern=r'\d{4}'):        """Internal helper to load temporal TIF stacks into a 3D numpy array (T, H, W)."""        files = sorted(            [f for f in os.listdir(var_path) if f.endswith(".tif")],             key=lambda x: int(re.findall(r"\d+", re.search(pattern, x).group())[0])        )        stack = []        for f in files:            with rasterio.open(os.path.join(var_path, f)) as src:                arr = src.read(1).astype(float)                if src.nodata is not None:                    arr[arr == src.nodata] = np.nan                stack.append(arr)                if self.profile is None:                    # Save profile from the first file to use for output GeoTIFFs                    self.profile = src.profile.copy()        return np.stack(stack)    def _calculate_trend(self, data_1d):        """Calculates the linear slope (trend) of a 1D time series."""        y = data_1d.flatten()        x = np.arange(len(y))        mask = ~np.isnan(y)        if np.sum(mask) < 2: return 0        slope, _, _, _, _ = stats.linregress(x[mask], y[mask])        return slope    def _normalize(self, data_1d):        """Min-Max Normalization: Xm = (x - min) / (max - min)."""        d_min, d_max = np.nanmin(data_1d), np.nanmax(data_1d)        if d_max == d_min: return np.zeros_like(data_1d)        return (data_1d - d_min) / (d_max - d_min)    def RR(self, y_var_name, pattern=r'\d{4}', k_best=2):        """        Main Ridge Regression function.        :param y_var_name: Folder name of the dependent variable (e.g., 'ET').        :param k_best: Regularization parameter (Lambda).        :return: Dictionary containing 'summary' (DataFrame) and 'maps' (2D Arrays).        """        ......    def save_results_to_tif(self, results, output_dir):        """Exports all result maps to GeoTIFF and summary to CSV."""        if not os.path.exists(output_dir):            os.makedirs(output_dir)        # Update metadata profile for output        out_profile = self.profile.copy()        out_profile.update({            "driver": "GTiff",            "dtype": 'float32',            "count": 1,            "nodata": self.nodata})        # Save each map        for name, data in results['maps'].items():            out_path = os.path.join(output_dir, f"{name}.tif")            out_data = data.astype(np.float32)            out_data[np.isnan(out_data)] = self.nodata            with rasterio.open(out_path, "w", **out_profile) as dst:                dst.write(out_data, 1)        # Save summary table        results['summary'].to_csv(os.path.join(output_dir, "summary.csv"), index=False)        print(f"Results successfully saved to {output_dir}")

import numpy as npimport matplotlib.pyplot as pltimport geosis.rdr as rdr

data_dir = r'E:\公众号\岭回归\data\vars'

# 创建数据集：data_dir 是所有变量文件夹的上一级文件夹inst = rdr.RidgeTrend(data_dir)# y_var_name 指的是因变脸名称，pattern 指的是每个变量中# 变量名称和时间名称的关系（例如‘\d{4}’代表变量名+年份# k_best 是岭回归正则化参数（需通过交叉验证等方法确定最优值）RidgeRe = inst.RR(y_var_name = 'SOT', pattern=r'\d{4}', k_best=2)

RidgeRe.get('summary')

maps = RidgeRe.get('maps')# 可视化im = plt.imshow(maps['r2'], cmap='RdYlGn')plt.colorbar(im)

inst.save_results_to_tif(RidgeRe, output_dir=r'E:\公众号\岭回归\Result')