基于Python代码采用地理探测器建成区驱动力分析

# -*- coding: utf-8 -*-

import pandas as pd

import numpy as np

import matplotlib.pyplot as plt

import warnings

from geodetector import Detector  # 核心地理探测器库

from sklearn.preprocessing import KBinsDiscretizer  # 连续变量离散化

warnings.filterwarnings('ignore')

# == 第一步：数据准备与预处理 ==

def load_and_preprocess_data(data_path):

    """

    加载建成区+驱动因子数据，完成预处理（缺失值、异常值、空间匹配）

    :param data_path: 数据文件路径（Excel/Csv，行=研究单元，列=建成区面积+驱动因子）

    :return: 因变量（建成区面积）、自变量（驱动因子）、因子名称列表

    """

    # 1. 读取数据（示例结构：行=区县/网格，列=建成区面积、夜光值、GDP、人口密度、坡度等）

    df = pd.read_csv(data_path, encoding='utf-8')  # Excel用pd.read_excel

    # 2. 定义变量（根据你的研究修改！）

    # 因变量：建成区范围表征（如建成区面积、建成区占比）

    target_col = '建成区面积_km2'

    # 自变量：建成区驱动因子（地理探测器候选因子）

    factor_cols = [

        '夜光强度均值', 'GDP_万元', '人口密度_人/km2',

        '坡度_mean', '距市中心距离_km', '第二产业占比_%'

    ]

    # 3. 数据清洗

    # 删除缺失值

    df = df.dropna(subset=[target_col] + factor_cols)

    # 异常值处理（3倍标准差法）

    for col in [target_col] + factor_cols:

        mean = df[col].mean()

        std = df[col].std()

        df = df[(df[col] >= mean - 3*std) & (df[col] <= mean + 3*std)]

    # 4. 提取变量

    y = df[target_col].values  # 因变量（建成区面积）

    X = df[factor_cols].values  # 自变量（驱动因子）

    print(f"✅ 数据预处理完成：有效样本数={len(df)}，驱动因子数={len(factor_cols)}")

    return y, X, factor_cols, df

# == 第二步：连续因子离散化（地理探测器核心要求） ==

def discretize_factors(X, factor_cols, n_bins=5):

    """

    地理探测器要求变量离散化，对连续因子进行分箱（等距/分位数）

    :param X: 自变量矩阵

    :param factor_cols: 因子名称

    :param n_bins: 分箱数（建议5-8箱）

    :return: 离散化后的自变量矩阵、分箱后的因子数据框

    """

    # 初始化分箱器（分位数分箱，更适配地理数据分布）

    discretizer = KBinsDiscretizer(

        n_bins=n_bins, 

        encode='ordinal',  # 编码为整数

        strategy='quantile'  # 分位数分箱（避免极端值影响）

    )

    X_discrete = discretizer.fit_transform(X)

    # 转换为DataFrame便于查看

    X_discrete_df = pd.DataFrame(X_discrete, columns=factor_cols)

    # 打印分箱结果（论文可用）

    print("\n📌 因子离散化分箱结果：")

    for col in factor_cols:

        bins_count = X_discrete_df[col].value_counts().sort_index()

        print(f"{col}：分箱数={len(bins_count)}，各箱样本数={bins_count.values}")

    return X_discrete, X_discrete_df

# ==第三步：地理探测器分析（因子+交互探测） ==

def run_geodetector(y, X_discrete, factor_cols):

    """

    运行地理探测器：因子探测（q值）、交互探测（交互类型）

    """

    # 初始化地理探测器

    detector = Detector()

    # 传入数据（因变量+离散化自变量）

    detector.fit(y=y, X=X_discrete, X_names=factor_cols)

    # 1. 因子探测：计算各因子的q值和显著性（核心）

    factor_detector = detector.detect()

    # 整理因子探测结果

    factor_result = pd.DataFrame({

        '驱动因子': factor_cols,

        'q值': [factor_detector[col]['q'] for col in factor_cols],

        'p值': [factor_detector[col]['p'] for col in factor_cols],

        '显著性': ['显著' if p < 0.05 else '不显著' for p in [factor_detector[col]['p'] for col in factor_cols]]

    }).sort_values(by='q值', ascending=False)

    # 2. 交互探测：分析因子间的交互效应

    interaction_detector = detector.interaction_detect()

    # 整理交互探测结果（筛选核心交互对）

    interaction_result = []

    for i in range(len(factor_cols)):

        for j in range(i+1, len(factor_cols)):

            factor1 = factor_cols[i]

            factor2 = factor_cols[j]

           interaction_type = interaction_detector.loc[factor1, factor2]['interaction_type']

          q_interaction = interaction_detector.loc[factor1, factor2]['q']

     q1 = factor_result[factor_result['驱动因子']==factor1]['q值'].values[0]

      q2 = factor_result[factor_result['驱动因子']==factor2]['q值'].values[0]

            interaction_result.append({

                '因子1': factor1,

                '因子2': factor2,

                'q1': q1,

                'q2': q2,

                'q交互': q_interaction,

                '交互类型': interaction_type

            })

interaction_result_df=pd.DataFrame(interaction_result).sort_values(by='q交互', ascending=False)

    # 打印结果

    print("\n🎯 因子探测结果（按q值降序）：")

    print(factor_result)

    print("\n🎯 交互探测结果（前5对核心交互）：")

    print(interaction_result_df.head(5))

    return factor_result, interaction_result_df, detector

# == 第四步：结果可视化（论文级图表） ==

def visualize_geodetector(factor_result, interaction_result_df, save_path):

    """

    可视化地理探测器结果：因子q值柱状图、交互探测热力图

    """

    plt.rcParams['font.sans-serif'] = ['SimHei']  # 中文显示

    plt.rcParams['axes.unicode_minus'] = False

    # 1. 绘制因子q值柱状图（核心）

    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(16, 6), dpi=300)

    # 因子q值图

    factor_result_sorted = factor_result.sort_values(by='q值')