# 空间分析与统计 | 关联：12_ArcGIS_Pro.md 21_Python_GIS生态.md 28_项目案例集.md | 来源：旧第十六篇

# 第十六篇：GIS数据分析与空间统计

> 来源：ArcGIS Pro 空间统计工具箱官方文档 + GeoDa 文档 + 豆包GIS专家知识框架
> 最后更新：2026-05-23 | 核心要点：理解空间非平稳性与空间自相关

---

## 16.1 空间统计核心概念

### 16.1.1 为什么需要空间统计？

传统统计学假设观测值相互独立，但地理数据（Tobler 第一定律）天然存在空间依赖：**"Everything is related to everything else, but near things are more related than distant things."**

| 对比维度 | 传统统计 | 空间统计 |
|---------|---------|---------|
| 核心假设 | 观测独立 | 空间自相关存在 |
| 研究目标 | 变量间关系 | 变量关系 + 空间格局 |
| 关键指标 | p值、R² | Moran's I、Gi* Z-score |
| 典型方法 | 回归、方差分析 | GWR、空间聚类、克里金 |
| 输出形式 | 数值结果 | 地图可视化 + 统计量 |

### 16.1.2 ArcGIS Pro 空间统计工具箱结构

ArcGIS Pro 空间统计工具箱（Spatial Statistics Tools）含 **6 个工具集**：

| 工具集 | 英文名 | 典型工具 | 用途 |
|--------|--------|---------|------|
| **分析模式** | Analyzing Patterns | 空间自相关(Moran's I)、高低聚类(Getis-Ord General G)、增量空间自相关 | 判断数据是否存在空间聚集 |
| **评估敏感性** | Assessing Sensitivity | 权重矩阵构建、空间自相关敏感性 | 评估参数选择对结果的稳定性 |
| **聚类分布制图** | Mapping Clusters | 热点分析(Getis-Ord Gi*)、聚类与异常值分析(Anselin Local Moran's I)、分组分析 | 识别聚集的具体位置 |
| **度量地理分布** | Measuring Geographic Distributions | 中心要素、方向分布(标准差椭圆)、平均中心 | 描述数据空间分布特征 |
| **空间关系建模** | Modeling Spatial Relationships | GWR、OLS、探索性回归 | 建立空间回归模型 |
| **空间组件实用程序** | Space Time Utilities | 时空立方体创建、时空热点分析 | 时空一体化分析 |

> **神经连接**：本节工具 → 第十二篇 ArcPy 参考手册（直接调用 `arcpy.stats` 模块）。

---

## 16.2 空间自相关分析（Moran's I）

### 16.2.1 全局 Moran's I 原理

**定义**：度量整个研究区域内空间自相关的总体程度。

**公式**：
```
        N     Σᵢ Σⱼ wᵢⱼ (xᵢ - x̄)(xⱼ - x̄)
I =  ─────── · ──────────────────────────
      Σᵢ Σⱼ wᵢⱼ        Σᵢ (xᵢ - x̄)²
```

其中：N = 要素数量，wᵢⱼ = 空间权重矩阵，xᵢ = 要素 i 的属性值，x̄ = 均值。

| I 值范围 | 含义 | 空间格局 |
|---------|------|---------|
| +1（接近） | 强正相关 | 高值聚在一起，低值聚在一起（聚集） |
| 0（接近） | 无自相关 | 随机分布 |
| -1（接近） | 强负相关 | 高值与低值交替（离散/棋盘格） |

**统计检验**：H₀ = 空间随机分布。Z 得分和 p 值判定显著性：
- |Z| > 1.96（p < 0.05）：显著
- |Z| > 2.58（p < 0.01）：非常显著

### 16.2.2 局部 Moran's I（LISA）

**Anselin Local Moran's I**：识别空间聚集的具体位置，每个观测值输出一个 I 值。

**Moran 散点图四象限**：

```
        Wz（邻域标准化值）
            ↑
    LH (Low-High)    │    HH (High-High)
    "洼地"被高值包围 │  "热点"高值聚集
    ─────────────────┼─────────────────→ Z（标准化值）
    LL (Low-Low)     │    HL (High-Low)
    "冷点"低值聚集   │  "孤岛"高值孤立
```

| 象限 | 类型 | 地理含义 | 示例 |
|------|------|---------|------|
| **HH** | 高-高聚集 | 自身高值，邻域也高值 | 市中心高房价聚集区 |
| **LL** | 低-低聚集 | 自身低值，邻域也低值 | 偏远地区低房价带 |
| **HL** | 高-低异常 | 自身高值，邻域低值 | 农村区域突然的高收入飞地 |
| **LH** | 低-高异常 | 自身低值，邻域高值 | 富人区里的低收入单元 |

**ArcPy 调用**：

```python
import arcpy
arcpy.stats.AnselinLocalMoransI(
    Input_Feature_Class="居民点",
    Input_Field="房价",
    Output_Feature_Class="LISA_Results.shp",
    Conceptualization_of_Spatial_Relationships="CONTIGUITY_EDGES_CORNERS",
    Distance_Method="EUCLIDEAN",
    Number_of_Neighbors=None
)
# 输出字段：LMiIndex(Local Moran's I)、LMiZScore、LMiPValue、COType(HH/HL/LH/LL)
```

### 16.2.3 双变量 Moran's I（Bivariate）

用于分析**一个变量的值与另一个变量邻域值**之间的空间关联：

```python
# 房价 与 邻域教育水平 的空间关联
arcpy.stats.BivariateMoransI(
    Input_Feature_Class="小区",
    Input_Field="房价",
    Related_Feature_Class="学校",
    Related_Field="升学率"
)
```

### 16.2.4 空间权重矩阵构建

空间权重矩阵 wᵢⱼ 决定了"邻近"的定义，直接影响 Moran's I 结果：

| 权重类型 | 构建方式 | 适用场景 |
|---------|---------|---------|
| **固定距离** | 指定距离阈值 d，≤d 为邻域 | 面状数据大小均匀 |
| **K 近邻** | 最近的 K 个要素为邻域 | 点数据密度不均 |
| **面邻接(共边)** | CONTIGUITY_EDGES_ONLY | 行政区域分析 |
| **面邻接(共边共点)** | CONTIGUITY_EDGES_CORNERS | 面状数据（推荐默认） |
| **反距离** | wᵢⱼ = 1/dᵢⱼᵖ | 距离衰减效应强 |

> **关键选择**：使用 `增量空间自相关` 工具确定最佳距离阈值 — 取 Moran's I Z-score 峰值对应的距离。

---

## 16.3 热点分析（Getis-Ord Gi*）

### 16.3.1 原理与公式

**Getis-Ord Gi*** 识别统计显著的高值簇（热点）和低值簇（冷点）：

```
        Σⱼ wᵢⱼ xⱼ - x̄ Σⱼ wᵢⱼ
Gi* = ─────────────────────────────
              √[S · ...]
```

**核心输出**：每个要素的 **Z 得分** 和 **p 值**。

### 16.3.2 结果解读

| Z 得分 | 置信度 | 分类 | 解读 |
|--------|--------|------|------|
| > +2.58 | 99% | Gi_Bin = +3 | **热点**：显著高值空间聚集 |
| +1.96 ~ +2.58 | 95% | Gi_Bin = +2 | 次热点 |
| +1.65 ~ +1.96 | 90% | Gi_Bin = +1 | 弱热点 |
| -1.65 ~ +1.65 | — | Gi_Bin = 0 | 不显著（随机） |
| -1.96 ~ -1.65 | 90% | Gi_Bin = -1 | 弱冷点 |
| -2.58 ~ -1.96 | 95% | Gi_Bin = -2 | 次冷点 |
| < -2.58 | 99% | Gi_Bin = -3 | **冷点**：显著低值空间聚集 |

**与 Moran's I 对比**：
| 特征 | Global Moran's I | Getis-Ord Gi* |
|------|-----------------|---------------|
| 分析尺度 | 全局（一个值） | 局部（每要素一个值） |
| 输出类型 | 整体模式判断 | 热点/冷点地图 |
| 统计量 | I 值 + Z-score + p-value | Z-score + p-value |
| 并行工具 | Anselin Local Moran's I | — |

### 16.3.3 ArcPy 调用

```python
arcpy.stats.HotSpots(
    Input_Feature_Class="犯罪点",
    Input_Field="案件数",
    Output_Feature_Class="Crime_HotSpots.shp",
    Conceptualization_of_Spatial_Relationships="FIXED_DISTANCE_BAND",
    Distance_Method="EUCLIDEAN",
    Distance_Band_or_Threshold_Distance="500 Meters",
    Self_Potential_Field=None
)
# 输出字段：GiZScore、GiPValue、Gi_Bin(-3~+3)
```

---

## 16.4 聚类分析方法

### 16.4.1 方法对比

| 方法 | 类型 | 原理 | 需预设组数 | ArcGIS 工具 |
|------|------|------|-----------|------------|
| **K-Means** | 非空间 | 最小化组内距离平方和 | ✅ 需要 | 分组分析(Grouping Analysis) |
| **DBSCAN** | 空间 | 基于密度的空间聚类 | ❌ 自动 | 密度聚类(Density-based Clustering) |
| **层次聚类** | 非空间 | 构建聚类树 | ✅ 需要 | 分组分析(选项) |
| **时空立方体** | 时空 | 三维立方体+时间序列 | ❌ | 时空模式挖掘工具集 |

### 16.4.2 DBSCAN 密度聚类

**关键参数**：
- **ε（eps）**：邻域半径
- **MinPts**：最小点数形成簇
- **核心点 / 边界点 / 噪声点** 三分法

```python
arcpy.stats.DensityBasedClustering(
    Input_Point_Features="POI",
    Output_Features="POI_Clusters.shp",
    Cluster_Method="DBSCAN",
    Min_Features_per_Cluster=5,
    Search_Distance="100 Meters",
    Cluster_Sensitivity=50
)
# 输出字段：CLUSTER_ID（噪声点为-1）、COLOR_ID（渲染用）
```

### 16.4.3 分组分析（Grouping Analysis）

多属性聚类，支持 K-Means / 层次聚类：

```python
arcpy.stats.GroupingAnalysis(
    Input_Feature_Class="社区",
    Unique_ID_Field="社区ID",
    Output_Feature_Class="Community_Groups.shp",
    Number_of_Groups=5,
    Analysis_Fields=["人口密度", "绿地率", "商业POI数", "房价均价"],
    Spatial_Constraints="CONTIGUITY_EDGES_CORNERS"  # 面邻接约束
)
```

---

## 16.5 地理加权回归（GWR）

### 16.5.1 为什么需要 GWR？

**普通最小二乘（OLS）** 假设变量关系在整个研究区恒定（空间平稳性）：
```
y = β₀ + β₁x₁ + β₂x₂ + ... + ε    （β₁ 全域恒定）
```

**GWR** 允许系数在空间上变化（空间非平稳性）：
```
y(uᵢ,vᵢ) = β₀(uᵢ,vᵢ) + β₁(uᵢ,vᵢ)x₁ + β₂(uᵢ,vᵢ)x₂ + ...   （β 随位置变化）
```

**示例**：房价与地铁距离的关系 — 市中心系数可能为 +2000 元/米，郊区可能仅 +500 元/米。

### 16.5.2 GWR 关键参数

| 参数 | 含义 | 选择方法 |
|------|------|---------|
| **带宽（Bandwidth）** | 影响每个局部回归的邻域范围 | AICc 自动优化 / 手动指定 |
| **核函数（Kernel）** | 权重衰减方式 | FIXED（固定距离）/ ADAPTIVE（自适应K近邻） |
| **AICc** | 模型拟合优度 | 越小越好，与 OLS 对比确定 GWR 是否显著改善 |

### 16.5.3 OLS → GWR 工作流

```
1. OLS 回归 → 检查 6 项诊断
   ├ 残差是否正态？       （直方图）
   ├ 残差有无空间自相关？  （空间自相关工具）
   ├ 自变量共线性？        （VIF < 7.5）
   ├ 残差异方差？          （Kohenker BP 统计量）
   ├ 模型偏差？            （Jarque-Bera 统计量）
   └ R² 和调整 R²
   
2. 如果残差空间自相关显著 → 存在空间非平稳性 → 使用 GWR

3. GWR 输出对比：
   ├ GWR R² > OLS R²？    （解释力提升）
   ├ GWR AICc < OLS AICc？（拟合改善 ≥ 3 点）
   └ 局部 R² 地图：哪里模型效果好？
```

### 16.5.4 ArcPy 调用

```python
# Step 1: OLS
arcpy.stats.OrdinaryLeastSquares(
    Input_Feature_Class="小区",
    Unique_ID_Field="ID",
    Output_Feature_Class="OLS_Results.shp",
    Dependent_Variable="房价",
    Explanatory_Variables=["距地铁", "面积", "房龄"]
)

# Step 2: GWR（如果 OLS 残差有空间自相关）
arcpy.stats.GeographicallyWeightedRegression(
    Input_Feature_Class="小区",
    Dependent_Variable="房价",
    Explanatory_Variables=["距地铁", "面积"],
    Output_Feature_Class="GWR_Results.shp",
    Kernel_type="ADAPTIVE",
    Band_Method="AICc",
    Output_Predicted_Values="GWR_Predict.shp"
)
# 输出字段：Local_R2、Predicted、Residual、各系数(C1_距地铁, C2_面积)、Intercept
```

### 16.5.5 探索性回归（Exploratory Regression）

自动测试所有解释变量组合，输出满足 OLS 诊断的组合列表：

```python
arcpy.stats.ExploratoryRegression(
    Input_Feature_Class="小区",
    Dependent_Variable="房价",
    Candidate_Explanatory_Variables=["距地铁", "面积", "房龄", "朝向", "楼层", "容积率"],
    Weights_Matrix_File="小区_权重.swm",
    Output_Report_File="E:/Exploratory_Report.txt",
    Maximum_Number_of_Explanatory_Variables=4
)
```

---

## 16.6 空间插值方法对比

### 16.6.1 三大方法原理

| 方法 | 原理 | 优点 | 缺点 | 适用场景 |
|------|------|------|------|---------|
| **IDW**（反距离权重） | ẑ = Σ(zᵢ/dᵢᵖ)/Σ(1/dᵢᵖ) | 简单直观，计算快 | 不平滑，边界效应 | 采样点密集均匀 |
| **克里金（Kriging）** | 基于变异函数的最优无偏估计 | 提供预测误差，理论最优 | 需半变异函数建模，计算量大 | 任何场景（地统计标准方法） |
| **样条（Spline）** | 通过所有采样点的最小曲率曲面 | 平滑美观 | 可能产生不合理的极值 | 地形/气象连续表面 |

### 16.6.2 克里金变体

| 变体 | 适用条件 | ArcGIS 工具 |
|------|---------|------------|
| **普通克里金**（Ordinary） | 均值未知但恒定 | `Kriging` (type=ORDINARY) |
| **简单克里金**（Simple） | 均值已知且恒定 | `Kriging` (type=SIMPLE) |
| **泛克里金**（Universal） | 存在趋势（与坐标相关） | `Kriging` (type=UNIVERSAL) |
| **指示克里金**（Indicator） | 分类/超过阈值概率 | `IndicatorKriging` |
| **协同克里金**（Cokriging） | 利用辅助变量提高精度 | `Cokriging` |
| **经验贝叶斯克里金**（EBK） | 自动化建模，无需手动调参 ⭐ | `EmpiricalBayesianKriging` |

### 16.6.3 半变异函数（Semivariogram）

**核心公式**：
```
γ(h) = ½ · E[(z(s) - z(s+h))²]
```

**三参数**：

```
γ(h)
 ↑
 │  ┌──────────────────────── 基台（Sill）= 总方差
 │ ╱
 │╱   ↗ 块金（Nugget）= 微观变异 + 测量误差
 └──────────────────────→ h
      ↑
     变程（Range）= 空间自相关距离
```

| 参数 | 含义 | 影响 |
|------|------|------|
| **块金（Nugget）** | h=0 处的变异 | 大 = 测量误差大或微观变异大 |
| **基台（Sill）** | 总方差 | 决定插值的不确定性上限 |
| **变程（Range）** | 自相关消失距离 | 超过此距离的采样点不参与估计 |

**常用模型**：球状（Spherical）、指数（Exponential）、高斯（Gaussian）、线性（Linear）。

### 16.6.4 插值对比速查

```python
# IDW - 快速插值
arcpy.sa.Idw(
    in_point_features="高程点",
    z_field="ELEV",
    cell_size=5,
    power=2  # 距离衰减指数，默认2
)

# 普通克里金
arcpy.sa.Kriging(
    in_point_features="高程点",
    z_field="ELEV",
    kriging_model="Spherical",  # 球状模型
    cell_size=5
)

# 经验贝叶斯克里金 - 推荐首选 ⭐
arcpy.sa.EmpiricalBayesianKriging(
    in_point_features="高程点",
    z_field="ELEV",
    cell_size=5,
    transformation_type="NONE"  # 或 EMPIRICAL / LOGEMPIRICAL
)

# Spline 样条
arcpy.sa.Spline(
    in_point_features="高程点",
    z_field="ELEV",
    cell_size=5,
    spline_type="REGULARIZED"  # 或 TENSION
)
```

### 16.6.5 交叉验证（Cross Validation）

评估插值精度的核心方法 — 每次留一个点，用其余点预测，比较实测 vs 预测：

| 指标 | 含义 | 理想值 |
|------|------|--------|
| **ME**（平均误差） | 偏差方向 | 接近 0 |
| **RMSE**（均方根误差） | 预测精度 | 越小越好 |
| **MSE**（标准化均方根） | 标准误差评估 | 接近 1 |
| **ASE**（平均标准误差） | 与 RMSE 接近 | ≈ RMSE |

---

## 16.7 空间采样设计

### 16.7.1 采样策略对比

| 策略 | 方法 | 优点 | 缺点 |
|------|------|------|------|
| **简单随机** | 纯随机抽取 | 统计基础好 | 可能遗漏空间覆盖率 |
| **系统采样** | 规则网格/等间距 | 空间覆盖均匀 | 可能错过周期模式 |
| **分层随机** | 先分层后随机 | 保证各类型代表性 | 需先验分层知识 |
| **空间均衡** | 最小化样本间空间自相关 | GIS 领域推荐 ⭐ | 需要专业工具 |
| **自适应采样** | 根据初步结果动态增加 | 高效利用资源 | 实施复杂 |

### 16.7.2 ArcGIS 空间采样工具

```python
# 创建随机点
arcpy.management.CreateRandomPoints(
    out_path="E:/sampling",
    out_name="RandomSamples.shp",
    constraining_feature_class="研究区",
    number_of_points_or_field=50
)

# 创建渔网（系统采样网格）
arcpy.management.CreateFishnet(
    out_feature_class="Grid_500m.shp",
    origin_coord="120.0 30.0",
    y_axis_coord="120.0 30.1",
    cell_width=500, cell_height=500,
    number_rows=0, number_columns=0,
    template="研究区"
)

# 子集要素（分层抽样）
arcpy.management.SubsetFeatures(
    in_features="土壤样点",
    out_training_feature_class="Train.shp",
    out_test_feature_class="Test.shp",
    subset_size_units="PERCENTAGE_OF_INPUT",
    subset_size=70
)
```

---

## 16.8 空间统计工作流最佳实践

```
1. 数据探索
   ├ 直方图 → 检查分布形态（正态？偏态？）
   ├ QQ图 → 正态性检验
   ├ 半变异函数云 → 全局趋势+各向异性
   └ 趋势分析 → 移除全局趋势

2. 全局模式判断
   ├ 空间自相关(Moran's I) → 有聚集？
   └ 高低聚类(Getis-Ord General G) → 高值还是低值聚集？

3. 局部格局定位
   ├ 热点分析(Getis-Ord Gi*) → 热点/冷点地图
   └ 聚类与异常值(LISA) → HH/HL/LL/LH 分类

4. 回归建模
   ├ OLS → 检查6项诊断
   ├ 探索性回归 → 筛选最优变量组合
   └ GWR → 局部系数地图（若空间非平稳）

5. 插值成面
   ├ 半变异函数建模 → 选择模型（球状/指数/高斯）
   ├ EBK → 自动最优插值
   └ 交叉验证 → 评估精度

6. 结果呈现
   ├ 显著性地图（仅显示 p<0.05 的区域）
   ├ 局部R²地图（GWR模型拟合度）
   └ 预测标准误差地图（Kriging插值不确定性）
```

---

## 16.9 可执行代码集 ⭐ V4.1 新增

### 16.9.1 空间权重矩阵自动构建与优选

```python
import arcpy

# 方案A: 增量空间自相关 — 自动确定最佳距离阈值
arcpy.stats.IncrementalSpatialAutocorrelation(
    Input_Features="房价点",
    Input_Field="price",
    Number_of_Distance_Bands=30,
    Beginning_Distance="100 Meters",
    Distance_Increment="100 Meters",
    Output_Report_File="E:/reports/ISA_Report.pdf",
    Output_Table="E:/results/ISA_Results.dbf"
)
# 读取报告 → 取 Z-score 峰值对应距离作为最佳带宽

# 方案B: 生成空间权重矩阵文件 (.swm)
arcpy.stats.GenerateSpatialWeightsMatrix(
    Input_Feature_Class="房价点",
    Unique_ID_Field="FID",
    Output_Spatial_Weights_Matrix_File="E:/results/KNN8.swm",
    Conceptualization_of_Spatial_Relationships="K_NEAREST_NEIGHBORS",
    Number_of_Neighbors=8,
    Row_Standardization="ROW_STANDARDIZATION"
)

# 方案C: 对比多种权重方案
weight_configs = [
    ("KNN4", "K_NEAREST_NEIGHBORS", 4),
    ("KNN8", "K_NEAREST_NEIGHBORS", 8),
    ("KNN16", "K_NEAREST_NEIGHBORS", 16),
    ("CONTIG", "CONTIGUITY_EDGES_CORNERS", None),
]

for name, method, k in weight_configs:
    out_swm = f"E:/results/{name}.swm"
    kwargs = {
        "Input_Feature_Class": "房价点",
        "Unique_ID_Field": "FID",
        "Output_Spatial_Weights_Matrix_File": out_swm,
        "Conceptualization_of_Spatial_Relationships": method,
        "Row_Standardization": "ROW_STANDARDIZATION"
    }
    if k:
        kwargs["Number_of_Neighbors"] = k
    arcpy.stats.GenerateSpatialWeightsMatrix(**kwargs)
    
    # 使用该权重计算 Moran's I 并比较 Z-score
    result = arcpy.stats.SpatialAutocorrelation(
        Input_Feature_Class="房价点",
        Input_Field="price",
        Conceptualization_of_Spatial_Relationships="GET_SPATIAL_WEIGHTS_FROM_FILE",
        Weights_Matrix_File=out_swm
    )
    print(f"{name}: Z={result.getOutput(0):.2f}, P={result.getOutput(1):.4f}")
```

### 16.9.2 异常值自动检测与剔除

```python
import arcpy
import numpy as np

# 方法1: 箱线图法 (IQR) — 轻量级
fc = "E:/data/points.shp"
field = "value"

# 读取属性到 numpy
arr = arcpy.da.FeatureClassToNumPyArray(fc, [field, "OID@"])
vals = arr[field]
q1, q3 = np.percentile(vals, [25, 75])
iqr = q3 - q1
lower = q1 - 1.5 * iqr
upper = q3 + 1.5 * iqr

outlier_oids = [r[1] for r in arr if r[0] < lower or r[0] > upper]
print(f"检测到 {len(outlier_oids)} 个异常值 (阈值: {lower:.2f}~{upper:.2f})")

# 标记异常值（不删除，仅标记）
arcpy.management.AddField(fc, "IS_OUTLIER", "SHORT")
with arcpy.da.UpdateCursor(fc, ["OID@", "IS_OUTLIER"]) as cursor:
    for row in cursor:
        row[1] = 1 if row[0] in outlier_oids else 0
        cursor.updateRow(row)

# 方法2: Local Moran's I 异常值检测（空间异常）
arcpy.stats.AnselinLocalMoransI(
    Input_Feature_Class=fc,
    Input_Field=field,
    Output_Feature_Class="E:/results/LISA_Outliers.shp",
    Conceptualization_of_Spatial_Relationships="CONTIGUITY_EDGES_CORNERS"
)
# COType = "HL" 或 "LH" → 空间异常值
```

### 16.9.3 GWR 带宽自动优选与对比

```python
import arcpy

# 多带宽对比测试
bandwidths = [200, 500, 1000, 2000, "AICc"]  # AICc=自动优选
results = []

for bw in bandwidths:
    out_fc = f"E:/results/GWR_BW{bw}.shp"
    bw_method = "AICc" if bw == "AICc" else "BANDWIDTH_PARAMETER"
    
    arcpy.stats.GeographicallyWeightedRegression(
        Input_Feature_Class="房价点",
        Dependent_Variable="price",
        Explanatory_Variables=["dist_metro", "area", "age"],
        Output_Feature_Class=out_fc,
        Kernel_type="ADAPTIVE",
        Band_Method=bw_method,
        Bandwidth=bw if bw != "AICc" else None
    )
    
    # 读取 AICc (从输出要素类的元数据)
    results.append({"bandwidth": bw, "fc": out_fc})

# 推荐：选择 AICc 最小的方案（自动包含在 AICc 计算中）

# 快速对比 OLS vs GWR
arcpy.stats.OrdinaryLeastSquares(
    Input_Feature_Class="房价点", Unique_ID_Field="FID",
    Output_Feature_Class="E:/results/OLS.shp",
    Dependent_Variable="price",
    Explanatory_Variables=["dist_metro", "area"]
)
# 对比 OLS 和 GWR 的 AICc：
# 若 GWR_AICc < OLS_AICc - 3 → GWR 显著改善
```

### 16.9.4 PySAL + GeoPandas 空间自相关

```python
import geopandas as gpd
import libpysal
from esda.moran import Moran, Moran_Local
from splot.esda import moran_scatterplot, lisa_cluster

# 读取数据
gdf = gpd.read_file("E:/data/districts.shp")

# 构建空间权重矩阵 (Queen 邻接)
w = libpysal.weights.Queen.from_dataframe(gdf)
w.transform = 'r'  # 行标准化

# 全局 Moran's I
mi = Moran(gdf['gdp_per_capita'], w)
print(f"Moran's I: {mi.I:.4f}, Z: {mi.z_sim:.2f}, P: {mi.p_sim:.4f}")
print(f"是否显著聚集: {'是' if mi.p_sim < 0.05 else '否'}")  # p < 0.05

# 局部 LISA
lisa = Moran_Local(gdf['gdp_per_capita'], w)
gdf['LISA_quadrant'] = lisa.q  # 1=HH, 2=LH, 3=LL, 4=HL
gdf['LISA_pvalue'] = lisa.p_sim

# 筛选显著区域 (p < 0.05)
gdf['significant'] = gdf['LISA_pvalue'] < 0.05

# 导出结果
gdf[['district_name', 'gdp_per_capita', 'LISA_quadrant', 'LISA_pvalue', 'geometry']].to_file(
    "E:/results/LISA_Results.shp"
)
```

---

> **V4.1 新增**：§16.9 可执行代码集（空间权重优选/异常值检测/GWR带宽对比/PySAL空间自相关）


<!-- wm:坤图_GIS:V1.0 -->
