点密度度量——计数和核密度

点密度度量——计数和核密度#

汇总运算对于聚合数据非常有用，无论是分析总体趋势还是可视化数据集。与简单地查看或显示地图上的点、线和多边形相比，汇总允许有效地分析和交流数据。本章将探讨两个强调数据集中的汇总运算：矩形或六边形网格中的计数点，或多边形和核密度的计数点。

import geopandas as gpd
# import geoplot as gplt
import math
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import rasterio
from rasterio.transform import Affine
from scipy import stats
from shapely.geometry import Polygon, box
from sklearn.datasets import fetch_species_distributions
from sklearn.neighbors import KernelDensity

计数点在矩形或六边形网格或多边形#

为了概述网格，创建新的多边形层，由网格和覆盖在另一个特征。可以在网格的每个单元格中总结该特性的一个方面。多边形层通常由渔网(矩形单元格)组成，但使用六边形作为网格正变得越来越普遍。

让我们定义一个函数，它将创建一个由指定边长的矩形或六边形组成的网格。

def create_grid(feature, shape, side_length, proj):
    '''创建由矩形或六边形组成的网格，其指定的边长（side_length）覆盖了输入特征的范围。'''

    # 通过创建缓冲区，稍微扩展特征范围的最小值和最大值
    # 降低特征直接落在 cell 边界（两个 cells 之间）的可能性
    # 缓冲区是投影相关的（归于单位）
    feature = feature.buffer(20)

    # 获取缓冲输入特征的范围
    min_x, min_y, max_x, max_y = feature.total_bounds

    # 创建空列表来保存将组成网格的单个单元格
    cells_list = []

    # 如果指定，创建方格网格
    if shape in ["square", "rectangle", "box"]:

        # 改编自 https://james-brennan.github.io/posts/fast_gridding_geopandas/
        # 创建并遍历将定义具有指定边长的列位置的 x 值列表
        for x in np.arange(min_x - side_length, max_x + side_length, side_length):
            # 创建并遍历将定义具有指定边长的行位置的 y 值列表
            for y in np.arange(min_y - side_length, max_y + side_length, side_length):
                # 创建具有指定边长的方框并添加到列表中
                cells_list.append(box(x, y, x + side_length, y + side_length))


    # 否则，创建六边形网格
    elif shape == "hexagon":
        # 设置水平位移，它将定义具有指定边长（基于 normal hexagon）的列位置。
        x_step = 1.5 * side_length

        # 设置垂直位移，将定义具有指定边长的行位置（基于 normal hexagon）
        # 这是两个堆叠在一起的六边形的中心之间的距离（垂直）
        y_step = math.sqrt(3) * side_length

        # 求中点（边的中心点到中点的距离，基于 normal hexagon）
        apothem = (math.sqrt(3) * side_length / 2)

        # 设置列号
        column_number = 0

        # 创建并遍历将定义具有垂直位移的列位置的 x 值列表
        for x in np.arange(min_x, max_x + x_step, x_step):
            # 创建并遍历将定义具有水平位移的列位置的 y 值列表
            for y in np.arange(min_y, max_y + y_step, y_step):
                # 创建具有指定边长的六边形
                hexagon = [[x + math.cos(math.radians(angle)) * side_length, y + math.sin(math.radians(angle)) * side_length] for angle in range(0, 360, 60)]
                # 在列表中添加六边形
                cells_list.append(Polygon(hexagon))

            # 检查列数是否为偶数
            if column_number % 2 == 0:
                # 如果是偶数，则将 y 的最小值和最大值展开，
                # 以垂直取代下一行展开值，以避免遗漏任何靠近特征范围的特征
                min_y -= apothem
                max_y += apothem

            # E否则，为奇数
            else:
                # 将 y 的最小值和最大值恢复到原始值
                min_y += apothem
                max_y -= apothem

            # 将列数增加1
            column_number += 1

    # 否则，触发异常
    else:
        raise Exception("Specify a rectangle or hexagon as the grid shape.")

    # 从单元格列表创建网格
    grid = gpd.GeoDataFrame(cells_list, columns = ['geometry'], crs = proj)

    # 创建一个列，为每个网格分配一个编号
    grid["Grid_ID"] = np.arange(len(grid))
    return grid

将通过计算网格中每个单元格中的井点数量来说明这种方法。可以通过两种不同的方式来实现这种方法，它们都有优点和缺点。

首先，将为这两个示例设置一些全局参数。

# 设置网格中单元格的边长
# 这取决于选择的投影，因为长度是在投影中指定的单位
side_length = 5000

# 设置网格形状
shape = "hexagon"
# shape = "rectangle"

分组#

此方法涉及使用空间连接将每个点分配到它所在的单元格。每个 cell 内的所有点随后被分组并计数。

目录