地理空间数据分析 | 核密度估计 | 基于邻居的查询

简介

本实验展示了一个基于邻居的地理空间数据查询示例（特别是核密度估计），使用基于哈弗辛距离度量构建的球树——即经纬度点之间的距离。数据集由菲利普斯等人（2006 年）提供。该示例使用底图库绘制南美洲的海岸线和国家边界。

虚拟机使用提示

虚拟机启动完成后，点击左上角切换到“笔记本”标签，以访问 Jupyter Notebook 进行练习。

有时，你可能需要等待几秒钟让 Jupyter Notebook 完成加载。由于 Jupyter Notebook 的限制，操作验证无法自动化。

如果你在学习过程中遇到问题，随时向 Labby 提问。课程结束后提供反馈，我们将立即为你解决问题。

导入所需库

第一步是导入本实验所需的库。在本实验中，我们将使用 numpy、matplotlib、fetch_species_distributions 和 KernelDensity 库。

import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import fetch_species_distributions
from sklearn.neighbors import KernelDensity

加载数据

下一步是加载菲利普斯等人（2006 年）提供的数据。该数据集包含两个物种，我们将用它们来演示核密度估计。

data = fetch_species_distributions()
species_names = ["Bradypus Variegatus", "Microryzomys Minutus"]

准备数据

现在我们将为核密度估计准备数据。我们将从数据集中提取纬度和经度信息，并将它们转换为弧度。

Xtrain = np.vstack([data["train"]["dd lat"], data["train"]["dd long"]]).T
ytrain = np.array(
    [d.decode("ascii").startswith("micro") for d in data["train"]["species"]],
    dtype="int",
)
Xtrain *= np.pi / 180.0  ## Convert lat/long to radians

构建网格

现在我们将从批次对象构建地图网格。我们将使用construct_grids函数来实现这一点。

def construct_grids(batch):
    """从批次对象构建地图网格

    参数
    ----------
    batch : 批次对象
        :func:`fetch_species_distributions` 返回的对象

    返回
    -------
    (xgrid, ygrid) : 一维数组
        与 batch.coverages 中的值对应的网格
    """
    ## 角单元格的 x、y 坐标
    xmin = batch.x_left_lower_corner + batch.grid_size
    xmax = xmin + (batch.Nx * batch.grid_size)
    ymin = batch.y_left_lower_corner + batch.grid_size
    ymax = ymin + (batch.Ny * batch.grid_size)

    ## 网格单元格的 x 坐标
    xgrid = np.arange(xmin, xmax, batch.grid_size)
    ## 网格单元格的 y 坐标
    ygrid = np.arange(ymin, ymax, batch.grid_size)

    return (xgrid, ygrid)

## 调用函数并将结果存储在 xgrid 和 ygrid 中
xgrid, ygrid = construct_grids(data)

准备数据网格

我们将为等高线图设置数据网格。我们将使用construct_grids函数来实现这一点。

X, Y = np.meshgrid(xgrid[::5], ygrid[::5][::-1])
land_reference = data.coverages[6][::5, ::5]
land_mask = (land_reference > -9999).ravel()

xy = np.vstack([Y.ravel(), X.ravel()]).T
xy = xy[land_mask]
xy *= np.pi / 180.0

绘制南美洲地图

现在我们将绘制南美洲地图，并标注出每个物种的分布情况。

fig = plt.figure()
fig.subplots_adjust(left=0.05, right=0.95, wspace=0.05)

for i in range(2):
    plt.subplot(1, 2, i + 1)

    print(" - 在球坐标中计算核密度估计")
    kde = KernelDensity(
        bandwidth=0.04, metric="haversine", kernel="gaussian", algorithm="ball_tree"
    )
    kde.fit(Xtrain[ytrain == i])

    Z = np.full(land_mask.shape[0], -9999, dtype="int")
    Z[land_mask] = np.exp(kde.score_samples(xy))
    Z = Z.reshape(X.shape)

    levels = np.linspace(0, Z.max(), 25)
    plt.contourf(X, Y, Z, levels=levels, cmap=plt.cm.Reds)

    if basemap:
        print(" - 使用 Basemap 绘制海岸线")
        m = Basemap(
            projection="cyl",
            llcrnrlat=Y.min(),
            urcrnrlat=Y.max(),
            llcrnrlon=X.min(),
            urcrnrlon=X.max(),
            resolution="c",
        )
        m.drawcoastlines()
        m.drawcountries()
    else:
        print(" - 根据覆盖范围绘制海岸线")
        plt.contour(
            X, Y, land_reference, levels=[-9998], colors="k", linestyles="solid"
        )
        plt.xticks([])
        plt.yticks([])

    plt.title(species_names[i])

plt.show()

总结

在本实验中，我们学习了如何对地理空间数据进行核密度估计。我们使用了菲利普斯等人（2006 年）的数据集来演示这项技术。我们还学习了如何使用 basemap 库绘制南美洲地图以及每个物种的分布情况。

物种分布的核密度估计

简介