在生态学和地理信息系统中，对地理空间数据进行核密度估计是一种常见的分析方法。本文将介绍如何使用基于哈维森距离度量的布朗树对地理空间数据进行核密度估计。哈维森距离是一种用于计算地球上两点之间距离的度量方法，特别适用于处理经纬度数据。
本例中使用的数据集由Phillips等人在2006年提供。如果可用，示例将使用basemap库来绘制南美洲的海岸线和国界。需要注意的是，本示例并不对数据进行任何学习（有关基于此数据集属性进行分类的示例，请参阅物种分布建模）。它仅展示了在地理坐标中观察到的数据点的核密度估计。
本例涉及的两个物种包括：
在进行核密度估计之前，首先需要构建地图网格。这可以通过从数据集中获取的Batch对象来完成。使用NumPy库中的arange函数来生成网格单元的x和y坐标。然后，使用这些坐标来构建一个网格，该网格对应于Batch对象中的覆盖值。
接下来，获取物种ID和位置的矩阵/数组。使用sklearn库中的fetch_species_distributions函数来获取数据，并从中提取物种名称。然后，将纬度和经度数据转换为弧度，以便于后续的计算。
为了在地图上绘制每个物种的分布，首先设置数据网格，以便为等高线图进行评估。使用NumPy库中的meshgrid函数来生成x和y的网格。然后，使用这些网格来评估核密度估计，并且只对陆地进行评估（-9999表示海洋）。
在绘制等高线图时，使用matplotlib库中的contourf函数来绘制密度等高线。如果basemap库可用，使用它来绘制南美洲的海岸线和国界。否则，使用数据集中的覆盖信息来绘制海岸线。
最后，使用matplotlib库中的show函数来显示地图。整个脚本的运行时间大约为3.409秒。
以下是一些与本示例相关的其他示例：
import matplotlib.pyplot as plt
import numpy as np
from sklearn.datasets import fetch_species_distributions
from sklearn.neighbors import KernelDensity
try:
from mpl_toolkits.basemap import Basemap
basemap = True
except ImportError:
basemap = False
def construct_grids(batch):
# 构建地图网格
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)
xgrid = np.arange(xmin, xmax, batch.grid_size)
ygrid = np.arange(ymin, ymax, batch.grid_size)
return xgrid, ygrid
# 获取物种ID和位置数据
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
# 构建等高线图的数据网格
xgrid, ygrid = construct_grids(data)
X, Y = np.meshgrid(xgrid[::5], ygrid[::5][:, ::-1])
land_reference = data.coverages[6][::5, ::5]
land_mask = (land_reference > -9999).ravel()
# 绘制南美洲地图及物种分布
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")
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()