Зональная статистика - Геоаналитика на Python

Зональная статистика — это метод пространственного анализа, позволяющий агрегировать значения растра внутри заданных векторных зон.

Иными словами, для каждой зоны (например, ячейки сетки или административной границы) рассчитываются статистические показатели значений пикселей, попадающих внутрь неё.

Этот метод позволяет перейти от непрерывной растровой поверхности (например, распределения населения) к агрегированным показателям по пространственным единицам.

В этом разделе мы применим зональную статистику к данным WorldPop и рассчитаем показатели по ячейкам регулярной сетки.

0. Импорт библиотек и подготовка данных¶

0.1 Импорт библиотек¶

import rasterio
import rasterstats as rs

import osmnx as ox
import geopandas as gpd
from shapely.geometry import box

0.2 Подготовка данных¶

Загрузим границу Тулы из OpenStreetMap

area_name = "Тула, Россия"

gdf_borders = ox.geocode_to_gdf(area_name)

gdf_borders.explore(tiles='cartodbpositron')

В предыдущем разделе мы подготовили растровый слой: обрезали его по границам города и перепроецировали в систему координат с метрическими единицами. Теперь загрузим полученный файл — tula_cropped_population_utm.tif — для дальнейшего анализа.

raster_path = "./data/tula_cropped_population_utm.tif"

raster = rasterio.open(raster_path)

1. Подготовка зон для анализа¶

В качестве зон будем использовать регулярную сетку, покрывающую территорию города Тулы.

1.1 Создание регулярной сетки¶

Сначала создадим функцию для генерации сетки заданного размера:

def create_regular_grid(gdf, square_size):
    # Определяем подходящую проекцию (в метрах)
    utm_crs = gdf.estimate_utm_crs()
    gdf_utm = gdf.to_crs(utm_crs)

    # Границы области
    minX, minY, maxX, maxY = gdf_utm.total_bounds

    grid_cells = []
    x, y = minX, minY

    # Генерация сетки
    while y <= maxY:
        while x <= maxX:
            grid_cells.append(box(x, y, x + square_size, y + square_size))
            x += square_size
        x = minX
        y += square_size

    fishnet = gpd.GeoDataFrame(geometry=grid_cells, crs=utm_crs)
    return fishnet

1.2 Построение сетки¶

Создадим сетку с шагом 1000 метров:

fishnet = create_regular_grid(gdf_borders, square_size=1000)
print(f"Создано ячеек: {len(fishnet)}")

Создано ячеек: 459

1.3 Визуализация сетки¶

fishnet.explore(tiles='cartodbpositron')

2. Подготовка данных для анализа¶

Перед расчётом зональной статистики важно убедиться, что растр и векторные данные находятся в одной системе координат.

2.1 Проверка систем координат¶

print("Grid CRS:", fishnet.crs)
print("Raster CRS:", raster.crs)

Grid CRS: EPSG:32637
Raster CRS: EPSG:32637

2.2 Приведение к одной системе координат¶

Если системы координат различаются, необходимо перепроецировать данные:

target_crs = raster.crs
# fishnet = fishnet.to_crs(target_crs)

3. Расчёт зональной статистики¶

Теперь рассчитаем статистику значений растра в пределах каждой ячейки сетки.

3.1 Расчёт суммы значений¶

Используем библиотеку rasterstats:

stats = rs.zonal_stats(
    fishnet,
    raster_path,
    stats="sum",
    geojson_out=True
)

3.2 Преобразование результата¶

gdf_stats = gpd.GeoDataFrame.from_features(stats, crs=raster.crs)
gdf_stats.head()

Теперь у нас есть GeoDataFrame, где для каждой ячейки рассчитана сумма значений растра.

4. Расчёт плотности населения¶

Рассчитаем плотность населения в каждой ячейке.

gdf_stats['density'] = gdf_stats['sum'] / (gdf_stats.geometry.area / 1_000_000)

Здесь:

sum — численность населения в ячейке
area — площадь ячейки (в м²)
делим на 1 000 000, чтобы получить км²

5. Визуализация результата¶

gdf_stats.explore(column='density', tiles='cartodbpositron')

Итог¶

В этом разделе мы познакомились с методом зональной статистики и применили его к данным WorldPop.

Мы:

создали регулярную сетку для исследуемой территории;
подготовили данные, приведя их к одной системе координат;
рассчитали сумму значений растра в пределах каждой ячейки;
вычислили плотность населения;
визуализировали полученный результат.

В результате мы перешли от непрерывного растрового слоя к агрегированным показателям, удобным для анализа и интерпретации.