Geospatial Data Science （5）: Spatial autocorrelation

1. 空间自相关和探索性空间数据分析

空间自相关与数据集中观察值的相似性与这些观察值的位置的相似性的程度有关。与传统的两个变量之间的相关关系（告知我们一个变量的值如何作为另一个变量的函数而变化）以及与它的时间序列对应物（将一个变量在某一特定时间点的值与以前的值联系起来）不完全一样，空间自相关将感兴趣的变量在一个特定地点的值与同一变量在周围地点的值联系起来。

这方面的一个关键概念是空间随机性：在这种情况下，观察的位置没有提供关于其价值的任何信息。换句话说，如果一个变量在空间上的分布没有明显的模式，那么它就是空间随机的。因此，空间自相关可以被正式定义为 “没有空间随机性”，这就为两类主要的自相关提供了空间，类似于传统的情况。正空间自相关，当类似的值倾向于聚集在类似的位置时；和负空间自相关，在类似的值倾向于分散和彼此之间的距离。

在本节课中，我们将学习如何在一个给定的数据集中探索空间自相关，询问数据的存在、性质和强度。为了做到这一点，我们将使用一套工具，统称为探索性空间数据分析（ESDA），专门为此目的而设计。ESDA方法的范围非常广泛，从不太复杂的方法，如折线图和一般的表格查询，到更先进和强大的方法，包括统计推理和对数据的地理维度的明确认识。本次会议的目的是让我们涉足后一组。

ESDA技术通常分为两大类：分析全球和本地空间自相关的工具。前者考虑数值位置所遵循的总体趋势，并对数据集中的*群集程度做出可能的陈述。数值在地理分布上是否普遍遵循一种特定的模式？*类似的数值与其他类似的数值是否比我们从纯粹的机会中预期的更接近？我们将通过使用Moran’s I统计学来练习全局空间自相关。

局部*空间自相关的工具则侧重于空间不稳定性：地图的部分区域偏离总体趋势。这里的意思是，即使在空间关联的性质和强度方面，数据有一个特定的趋势，但一些特定的地区可以与一般模式有相当大的差异。不管数值的总体集中程度如何，我们可以观察到靠近其他高（低）值的小块异常高（低）值，我们称之为热（冷）点。此外，也有可能观察到一些被低（高）值包围的高（低）值，我们将命名这些为 “空间异常值”。本节课我们将回顾探讨局部空间自相关的主要技术是局部空间关联指标（LISA）。

import seaborn as sns
import pandas as pd
import esda
from pysal.lib import weights
from splot.esda import (moran_scatterplot, lisa_cluster, plot_local_autocorrelation, plot_moran
)
from splot.libpysal import plot_spatial_weights
import geopandas as gpd
import numpy as np
import contextily as ctx
import matplotlib.pyplot as plt

2. Data

在这次会议上，我们将使用2016年公投中离开欧盟的投票结果，在地方当局层面。特别是，我们将专注于投票离开的空间分布，最终获胜。从技术角度看，你将与多边形打交道，这些多边形上有一个数值（投票给脱欧的选民的百分比）。

为方便起见，所有必要的数据都被汇集到一个文件中，其中包含英格兰、威尔士和苏格兰每个地方当局的地理信息，以及投票属性。该文件采用地理空间格式GeoPackage，与更传统的shapefile相比有几个优点（其中主要是需要一个文件而不是几个）。该文件可从课程网站下载。

# Read the file in
br = gpd.read_file("files/brexit.gpkg")

br.crs

<Geographic 2D CRS: EPSG:4326>
Name: WGS 84
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: World
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

现在，让我们在地方当局的ID上建立索引，同时将这些ID也保留为一列。

# Index table on the LAD ID
br = br.set_index("lad16cd", drop=False)
# Display summary
br.info()

<class 'geopandas.geodataframe.GeoDataFrame'>
Index: 380 entries, E06000001 to W06000024
Data columns (total 5 columns):#   Column     Non-Null Count  Dtype
---  ------     --------------  -----   0   objectid   380 non-null    int64   1   lad16cd    380 non-null    object  2   lad16nm    380 non-null    object  3   Pct_Leave  380 non-null    float64 4   geometry   380 non-null    geometry
dtypes: float64(1), geometry(1), int64(1), object(2)
memory usage: 17.8+ KB

3. 准备数据

# Plot polygons
ax = br.plot(alpha=0.5, color='red');
# Add background map, expressing target CRS so the basemap can be reprojected
ctx.add_basemap(ax, crs=br.crs)

3.1 空间权重矩阵

正如前面所讨论的，空间权重矩阵是将地理空间正式编码为数字形式的方式，因此它很容易被计算机（或统计方法）理解。我们已经看到了许多可以定义空间权重矩阵的概念性方式，如毗连性、基于距离或区块。

在这个例子中，我们将展示如何建立一个queen毗连矩阵，它认为如果两个观测值至少有一个共同的边界点，那么它们就是邻居。换句话说，数据集中的一对地方当局要想在这个 W W W下被认为是邻居，它们需要共享边界，或者换句话说，在某种程度上接触对方。

从技术上讲，我们将以与上次相同的方式建立毗连矩阵。我们将从GeoDataFrame开始，并将其传递给 PySAL中的女王毗连权重生成器（ps.weights.Queen.from_dataframe）。我们还将确保我们的数据表之前是以地方当局代码为索引的，因此 W W W也是以该表格为索引的。

# Create the spatial weights matrix
%time w = weights.Queen.from_dataframe(br, idVariable="lad16cd")

现在，我们拥有的w对象与我们过去创建的任何其他对象属于同一类型。因此，我们可以以相同的方式检查它。例如，我们可以检查谁是观察点E08000012的邻居：

w['E08000012']

{'E08000011': 1.0, 'E08000014': 1.0, 'E06000006': 1.0}

然而，我们计算 W W W的单元格返回了一个关于 "岛屿 "的警告。请记住，这些岛屿不一定是地理意义上的（尽管其中一些会是），而是数学意义上的：不与任何其他国家共享边界的地方当局，因此没有任何邻居。我们可以对它们进行检查并绘制地图，以更好地了解我们正在处理的问题。

ax = br.plot(color='k', figsize=(9, 9))
br.loc[w.islands, :].plot(color='red', ax=ax)