Алгоритм для нахождения минимального набора почтовых индексов для почти полного покрытия страны

Этот пост показывает, как выбрать минимальный набор почтовых индексов (PLZ) для почти полного покрытия страны, используя Python и некоторые распространённые библиотеки. Пример использует Германию, но подход можно адаптировать для других стран с аналогичными системами почтовых индексов.

Модель, которую мы здесь используем, заключается в том, что вокруг каждого выбранного почтового индекса (= центра) рисуется круг постоянного радиуса, и цель — покрыть как можно большую часть площади страны как можно меньшим количеством кругов. Это классическая задача оптимизации, часто называемая “задачей о покрытии множества”, которая является NP-трудной. Поэтому мы используем жадный аппроксимационный алгоритм с некоторыми улучшениями для нахождения хорошего решения за разумное время.

Этот алгоритм можно легко адаптировать для других стран, а также для других наборов POI.

Это в первую очередь стохастический алгоритм, поэтому он может не дать абсолютно оптимального результата, но должен быть очень близок к оптимальному в большинстве случаев. Кроме того, он выполняется за секунды даже на стандартном ноутбуке (с параметрами, как в примере ниже).

Поскольку этот код считывает таблицу PLZ как файл parquet, см. наш предыдущий пост Postleitzahlen und Koordinaten von GeoNames parsen для кода и запустите его с python parse_de.py -p postleitzahlen.parquet для создания необходимого файла.

Результаты

Используя пример ниже, мы в итоге выбираем 429 центров (почтовых индексов) для покрытия 99.5% Германии с радиусом 35 км вокруг каждого выбранного почтового индекса.

PLZ Coverage Germany.avif

Полный исходный код

plz-coverage.py

__author__ = "Uli Köhler"
__copyright__ = "Copyright 2025, Uli Köhler"
__license__ = "CC0 1.0 Universal (CC0 1.0) Public Domain Dedication"

import geopandas as gpd
import pandas as pd
import numpy as np
from shapely.geometry import Point
from shapely.ops import unary_union
import matplotlib.pyplot as plt
from tqdm.auto import tqdm
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.feature import ShapelyFeature
from sklearn.neighbors import KDTree
import cartopy.io.shapereader as shpreader

# Встроенные загрузчики PLZ и страны (самодостаточные)
def load_plz_data(plz_file_path: str) -> pd.DataFrame:
    """Загрузка табличных данных почтовых индексов из пути к файлу (parquet или CSV)."""
    if plz_file_path.endswith('.parquet') or plz_file_path.endswith('.parq'):
        return pd.read_parquet(plz_file_path)
    if plz_file_path.endswith('.csv') or plz_file_path.endswith('.txt'):
        return pd.read_csv(plz_file_path)
    # Попытка использования распространённых читателей как запасной вариант
    try:
        return pd.read_parquet(plz_file_path)
    except Exception:
        try:
            return pd.read_csv(plz_file_path)
        except Exception as exc:
            raise RuntimeError(f'Не удалось прочитать данные PLZ из {plz_file_path}: {exc}') from exc


def create_plz_points(plz_df: pd.DataFrame) -> gpd.GeoDataFrame:
    """Создание GeoDataFrame точек PLZ из DataFrame.
    Ожидается долгота/широта в столбцах с именами из ('lon','lng','longitude') и ('lat','latitude').
    Также гарантирует наличие столбца 'PLZ' путём копирования распространённых столбцов почтовых индексов или использования индекса как запасного варианта.
    """
    df = plz_df.copy()
    # толерантный поиск столбцов для координат
    lon_cols = [c for c in df.columns if c.lower() in ('lon', 'lng', 'longitude')]
    lat_cols = [c for c in df.columns if c.lower() in ('lat', 'latitude')]
    if not lon_cols or not lat_cols:
        raise ValueError('PLZ dataframe must contain lon and lat columns (e.g. lon, lat)')
    lon_col = lon_cols[0]
    lat_col = lat_cols[0]
    df = df.dropna(subset=[lon_col, lat_col]).copy()
    geometry = [Point(xy) for xy in zip(df[lon_col].astype(float), df[lat_col].astype(float))]
    gdf = gpd.GeoDataFrame(df, geometry=geometry, crs="EPSG:4326")
    # Сохранение стандартных имён столбцов для последующего кода
    if 'lon' not in gdf.columns:
        gdf['lon'] = gdf[lon_col].astype(float)
    if 'lat' not in gdf.columns:
        gdf['lat'] = gdf[lat_col].astype(float)

    # Гарантия наличия столбца PLZ (почтовый индекс). Сначала пробуем распространённые имена, иначе используем индекс как строку
    plz_cols = [c for c in gdf.columns if c.lower() in ('plz', 'postal_code', 'postal', 'postleitzahl')]
    if plz_cols:
        gdf['PLZ'] = gdf[plz_cols[0]].astype(str)
    else:
        # Если есть числовой столбец с именем вроде 'zip' или 'postcode', выбираем его
        zip_cols = [c for c in gdf.columns if c.lower() in ('zip', 'postcode')]
        if zip_cols:
            gdf['PLZ'] = gdf[zip_cols[0]].astype(str)
        else:
            # запасной вариант: создание PLZ из индекса для стабильного идентификатора
            gdf = gdf.reset_index(drop=False)
            gdf['PLZ'] = gdf['index'].astype(str)
            gdf.drop(columns=['index'], inplace=True)

    return gdf


def load_germany_from_natural_earth() -> object:
    """Чтение 10m Natural Earth стран и возврат геометрии Германии (shapely)"""
    # Использование пути shapereader cartopy и geopandas для загрузки 10m admin_0_countries shapefile
    shp_path = shpreader.natural_earth(resolution='10m', category='cultural', name='admin_0_countries')
    countries = gpd.read_file(shp_path)
    # пробуем распространённые имена ISO-столбцов для сопоставления Германии
    iso_cols = [c for c in ('ISO_A3', 'ISO3', 'adm0_a3', 'ADM0_A3', 'iso_a3') if c in countries.columns]
    if iso_cols:
        germany = countries[countries[iso_cols[0]] == 'DEU']
    else:
        # запасной вариант: пробуем сопоставление по NAME или SOVEREIGNT
        name_cols = [c for c in ('NAME_EN', 'NAME', 'SOVEREIGNT', 'ADMIN') if c in countries.columns]
        if name_cols:
            germany = countries[countries[name_cols[0]].str.contains('Germany', case=False, na=False)]
        else:
            raise RuntimeError('Не удалось определить ISO или столбец имени в странах Natural Earth')
    if germany.empty:
        raise RuntimeError('Геометрия Германии не найдена в странах Natural Earth')
    return germany.iloc[0].geometry


def filter_plz_in_germany(plz_gdf, germany_geometry):
    """Фильтрация точек PLZ, находящихся в пределах Германии"""
    plz_gdf = plz_gdf.to_crs("EPSG:3857")  # Web Mercator для вычисления расстояний
    germany_geometry = gpd.GeoSeries([germany_geometry], crs="EPSG:4326").to_crs("EPSG:3857")[0]

    # Создание буфера вокруг Германии для включения PLZ около границы
    buffered_germany = germany_geometry.buffer(50000)  # буфер 50 км

    # Фильтрация PLZ в пределах буферизованной Германии
    plz_in_germany = plz_gdf[plz_gdf.intersects(buffered_germany)]
    return plz_in_germany


def calculate_coverage(selected_plz, germany_geometry, radius_km=30):
    """Вычисление процента Германии, покрытого выбранными PLZ"""
    radius_m = radius_km * 1000
    coverage = unary_union(selected_plz.geometry.buffer(radius_m))
    germany_geometry = gpd.GeoSeries([germany_geometry], crs="EPSG:4326").to_crs("EPSG:3857")[0]
    covered_area = germany_geometry.intersection(coverage).area
    total_area = germany_geometry.area
    return covered_area / total_area * 100


def _generate_grid_points_in_polygon(polygon_3857, step_m: float) -> np.ndarray:
    """Генерация регулярной сетки точек (x, y) внутри данного полигона (EPSG:3857)."""
    xmin, ymin, xmax, ymax = polygon_3857.bounds
    xs = np.arange(xmin, xmax + step_m, step_m)
    ys = np.arange(ymin, ymax + step_m, step_m)
    points = []
    # Простой цикл подходит для ~20k-60k точек
    for x in xs:
        for y in ys:
            p = Point(x, y)
            if polygon_3857.contains(p):
                points.append((x, y))
    if not points:
        return np.empty((0, 2), dtype=float)
    return np.array(points, dtype=float)


def select_plz_max_coverage(
    plz_gdf: gpd.GeoDataFrame,
    germany_geometry,
    radius_km: float = 30,
    coverage_target: float = 0.99,
    sample_size: int = 300,
    grid_step_km: float = 5.0,
    exclusion_factor: float = 0.75,
    max_steps: int = 10000,
    patience: int = 500,
    validate_every: int = 50,
    random_state: int = 42,
    verbose: bool = True
) -> gpd.GeoDataFrame:
    """
    Единый алгоритм выбора точек PLZ с максимизацией покрытой площади (приближённо).

    Использует регулярную сетку точек внутри Германии (EPSG:3857) как proxy площади.
    Жадный рандомизированный: на каждом шаге выбирается подмножество оставшихся кандидатов
    и выбирается тот, который покрывает больше всего непокрытых точек сетки в радиусе.

    Улучшения:
    - Периодические точные проверки покрытия (unary_union) для избежания преждевременной остановки.
    - Радиус исключения для подавления кластерных выборов: после выбора центра,
      отбрасываются все оставшиеся PLZ в пределах (exclusion_factor * radius_km).

    Возвращает: GeoDataFrame (EPSG:3857) выбранных точек PLZ.
    """
    rng = np.random.default_rng(random_state)

    # Гарантия проецированной CRS для расстояния/площади и подготовка геометрии
    work_gdf = plz_gdf.to_crs("EPSG:3857").copy()
    germany_3857 = gpd.GeoSeries([germany_geometry], crs="EPSG:4326").to_crs("EPSG:3857")[0]

    radius_m = float(radius_km) * 1000.0
    exclusion_m = float(exclusion_factor) * radius_m
    grid_step_m = float(grid_step_km) * 1000.0

    # 1) Построение точек сетки внутри Германии как proxy покрытия
    if verbose:
        print(f"Построение сетки {grid_step_km:.1f} км над Германией для аппроксимации покрытия...")
    grid_xy = _generate_grid_points_in_polygon(germany_3857, step_m=grid_step_m)
    if grid_xy.shape[0] == 0:
        raise RuntimeError("Не удалось сгенерировать точки сетки внутри полигона Германии.")
    if verbose:
        print(f"Точки сетки: {grid_xy.shape[0]:,}")

    # 2) KDTree по сетке для быстрых запросов в радиусе
    tree_grid = KDTree(grid_xy, leaf_size=40, metric='euclidean')

    # 3) Предварительное вычисление индексов соседних точек сетки для каждого кандидата PLZ
    plz_xy = np.column_stack([work_gdf.geometry.x.values, work_gdf.geometry.y.values])
    if verbose:
        print("Предварительное вычисление окрестностей покрытия кандидат-сетка...")
    neighbor_lists = tree_grid.query_radius(plz_xy, r=radius_m, return_distance=False)

    # 3b) KDTree по кандидатам PLZ для фильтрации исключений
    tree_plz = KDTree(plz_xy, leaf_size=40, metric='euclidean')

    # 4) Жадный рандомизированный цикл выбора по маске непокрытой сетки
    uncovered = np.ones(grid_xy.shape[0], dtype=bool)
    selected_mask = np.zeros(len(work_gdf), dtype=bool)
    active_mask = np.ones(len(work_gdf), dtype=bool)  # кандидаты всё ещё доступные
    no_improve_steps = 0

    pbar = tqdm(total=max_steps, disable=not verbose, desc="Выбор центров (макс-покрытие)")
    for step in range(max_steps):
        pbar.update(1)

        # Периодическая точная проверка покрытия
        if step % validate_every == 0 and selected_mask.any():
            selected_tmp = work_gdf.loc[selected_mask]
            exact_cov = calculate_coverage(selected_tmp, germany_geometry, radius_km) / 100.0
            if step % (validate_every * 2) == 0:
                pbar.set_postfix_str(f"exact~{exact_cov*100:.1f}%")
            if exact_cov >= coverage_target:
                break

        # Приближённая доля покрытия по сетке
        covered_frac = 1.0 - (uncovered.sum() / uncovered.size)
        if covered_frac >= coverage_target and selected_mask.any():
            selected_tmp = work_gdf.loc[selected_mask]
            exact_cov = calculate_coverage(selected_tmp, germany_geometry, radius_km) / 100.0
            if exact_cov >= coverage_target:
                break

        remaining_idx = np.flatnonzero(active_mask)
        if remaining_idx.size == 0:
            break

        # Выборка кандидатов из активного множества
        k = int(min(sample_size, remaining_idx.size))
        sample = rng.choice(remaining_idx, size=k, replace=False) if k > 0 else []

        # Оценка предельного выигрыша по сетке
        best_gain = 0
        best_idx = None
        for idx in sample:
            neigh = neighbor_lists[idx]
            if neigh.size == 0:
                continue
            gain = int(uncovered[neigh].sum())
            if gain > best_gain:
                best_gain = gain
                best_idx = idx

        if best_idx is None or best_gain == 0:
            no_improve_steps += 1
            if no_improve_steps >= patience:
                break
            else:
                continue

        # Принятие лучшего кандидата
        selected_mask[best_idx] = True
        # Пометка точек сетки как покрытых
        uncovered[neighbor_lists[best_idx]] = False
        # Исключение близлежащих PLZ для уменьшения кластеризации
        if exclusion_m > 0:
            near = tree_plz.query_radius(plz_xy[[best_idx]], r=exclusion_m, return_distance=False)[0]
            active_mask[near] = False
        # Также удаление выбранного индекса
        active_mask[best_idx] = False
        no_improve_steps = 0

    pbar.close()

    selected = work_gdf.loc[selected_mask].copy()
    selected = selected.set_crs("EPSG:3857")
    return selected


def visualize_coverage(selected_plz, germany_geometry, radius_km=30):
    """Визуализация покрытия с использованием Cartopy"""
    # Конвертация в WGS84 для визуализации
    selected_plz = selected_plz.to_crs("EPSG:4326")
    germany_geometry = gpd.GeoSeries([germany_geometry], crs="EPSG:4326").to_crs("EPSG:4326")[0]

    # Создание фигуры с проекцией Cartopy
    fig = plt.figure(figsize=(12, 12))
    ax = fig.add_subplot(1, 1, 1, projection=ccrs.EqualEarth())

    # Добавление элементов карты
    ax.add_feature(cfeature.LAND, facecolor='lightgray')
    ax.add_feature(cfeature.OCEAN, facecolor='lightblue')
    ax.add_feature(cfeature.COASTLINE, linewidth=0.5)
    ax.add_feature(cfeature.BORDERS, linestyle=':', linewidth=0.5)

    # Создание ShapelyFeature для Германии
    germany_feature = ShapelyFeature([germany_geometry], ccrs.PlateCarree(), facecolor='none', edgecolor='red', linewidth=1)
    ax.add_feature(germany_feature)

    # Создание и отрисовка кругов (теперь ВСЕ PLZ)
    radius_m = radius_km * 1000
    buffers_3857 = gpd.GeoSeries(selected_plz.geometry).to_crs("EPSG:3857").buffer(radius_m)
    buffers_4326 = gpd.GeoSeries(buffers_3857, crs="EPSG:3857").to_crs("EPSG:4326")

    # Отрисовка каждого буфера по отдельности, чтобы перекрытия отображались как слоистые прозрачности
    for geom in buffers_4326:
        ax.add_geometries([geom], crs=ccrs.PlateCarree(), facecolor='blue', edgecolor='none', alpha=0.2, zorder=2)

    # Отрисовка выбранных точек PLZ поверх
    ax.scatter(selected_plz.geometry.x, selected_plz.geometry.y, color='red', s=5, transform=ccrs.PlateCarree(), label='Выбранные PLZ', zorder=3)

    # Установка охвата для Германии с отступами
    ax.set_extent([4, 16, 47, 56], crs=ccrs.PlateCarree())

    # Добавление линий сетки
    gl = ax.gridlines(draw_labels=True, linestyle='--')
    gl.top_labels = False
    gl.right_labels = False

    plt.title(f"Покрытие PLZ Германии с {len(selected_plz)} центрами (радиус {radius_km} км)")
    plt.legend(loc='upper right')
    plt.show()


def run_analysis(plz_file_path, radius_km=30, coverage_target=0.99):
    """Главная функция для анализа с использованием единого алгоритма макс-покрытия"""
    # Загрузка геометрии Германии из Natural Earth
    print("Загрузка геометрии Германии из Natural Earth...")
    germany_geometry = load_germany_from_natural_earth()

    # Загрузка и обработка данных PLZ
    print("Загрузка данных PLZ...")
    plz_df = load_plz_data(plz_file_path)
    print(f"Загружено {len(plz_df)} немецких почтовых индексов")

    plz_gdf = create_plz_points(plz_df)

    print("Фильтрация PLZ в пределах Германии...")
    plz_in_germany = filter_plz_in_germany(plz_gdf, germany_geometry)
    print(f"Найдено {len(plz_in_germany)} PLZ в пределах или около Германии")

    # Выбор центров PLZ с использованием единого алгоритма
    print("Выбор центров PLZ с алгоритмом макс-покрытия...")
    selected_plz = select_plz_max_coverage(
        plz_in_germany, germany_geometry, radius_km=radius_km, coverage_target=coverage_target, sample_size=300, grid_step_km=2.0, exclusion_factor=0.75, max_steps=10000, patience=500, validate_every=50, random_state=42, verbose=True
    )

    # Вычисление покрытия (точное, единое объединение)
    coverage_percent = calculate_coverage(selected_plz, germany_geometry, radius_km)
    print(f"\nВыбрано {len(selected_plz)} центров PLZ, покрывающих {coverage_percent:.1f}% Германии")

    # Визуализация покрытия
    visualize_coverage(selected_plz, germany_geometry=germany_geometry, radius_km=radius_km)

    # Защитный возврат результатов: включаются только существующие столбцы
    sel = selected_plz.to_crs("EPSG:4326")
    desired = ['PLZ', 'lat', 'lon', 'geometry']
    available = [c for c in desired if c in sel.columns or c == 'geometry']
    # Гарантия наличия геометрии
    if 'geometry' not in available:
        available.append('geometry')

    return {
        'selected_plz': sel[available],
        'coverage_percent': coverage_percent,
        'total_plz_used': len(selected_plz)
    }

__author__ = "Uli Köhler"
__copyright__ = "Copyright 2025, Uli Köhler"
__license__ = "CC0 1.0 Universal (CC0 1.0) Public Domain Dedication"

import geopandas as gpd
import pandas as pd
import numpy as np
from shapely.geometry import Point
from shapely.ops import unary_union
import matplotlib.pyplot as plt
from tqdm.auto import tqdm
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.feature import ShapelyFeature
from sklearn.neighbors import KDTree
import cartopy.io.shapereader as shpreader

# Встроенные загрузчики PLZ и страны (самодостаточные)
def load_plz_data(plz_file_path: str) -> pd.DataFrame:
    """Загрузка табличных данных почтовых индексов из пути к файлу (parquet или CSV)."""
    if plz_file_path.endswith('.parquet') or plz_file_path.endswith('.parq'):
        return pd.read_parquet(plz_file_path)
    if plz_file_path.endswith('.csv') or plz_file_path.endswith('.txt'):
        return pd.read_csv(plz_file_path)
    # Попытка использования распространённых читателей как запасной вариант
    try:
        return pd.read_parquet(plz_file_path)
    except Exception:
        try:
            return pd.read_csv(plz_file_path)
        except Exception as exc:
            raise RuntimeError(f'Не удалось прочитать данные PLZ из {plz_file_path}: {exc}') from exc


def create_plz_points(plz_df: pd.DataFrame) -> gpd.GeoDataFrame:
    """Создание GeoDataFrame точек PLZ из DataFrame.
    Ожидается долгота/широта в столбцах с именами из ('lon','lng','longitude') и ('lat','latitude').
    Также гарантирует наличие столбца 'PLZ' путём копирования распространённых столбцов почтовых индексов или использования индекса как запасного варианта.
    """
    df = plz_df.copy()
    # толерантный поиск столбцов для координат
    lon_cols = [c for c in df.columns if c.lower() in ('lon', 'lng', 'longitude')]
    lat_cols = [c for c in df.columns if c.lower() in ('lat', 'latitude')]
    if not lon_cols or not lat_cols:
        raise ValueError('PLZ dataframe must contain lon and lat columns (e.g. lon, lat)')
    lon_col = lon_cols[0]
    lat_col = lat_cols[0]
    df = df.dropna(subset=[lon_col, lat_col]).copy()
    geometry = [Point(xy) for xy in zip(df[lon_col].astype(float), df[lat_col].astype(float))]
    gdf = gpd.GeoDataFrame(df, geometry=geometry, crs="EPSG:4326")
    # Сохранение стандартных имён столбцов для последующего кода
    if 'lon' not in gdf.columns:
        gdf['lon'] = gdf[lon_col].astype(float)
    if 'lat' not in gdf.columns:
        gdf['lat'] = gdf[lat_col].astype(float)

    # Гарантия наличия столбца PLZ (почтовый индекс). Сначала пробуем распространённые имена, иначе используем индекс как строку
    plz_cols = [c for c in gdf.columns if c.lower() in ('plz', 'postal_code', 'postal', 'postleitzahl')]
    if plz_cols:
        gdf['PLZ'] = gdf[plz_cols[0]].astype(str)
    else:
        # Если есть числовой столбец с именем вроде 'zip' или 'postcode', выбираем его
        zip_cols = [c for c in gdf.columns if c.lower() in ('zip', 'postcode')]
        if zip_cols:
            gdf['PLZ'] = gdf[zip_cols[0]].astype(str)
        else:
            # запасной вариант: создание PLZ из индекса для стабильного идентификатора
            gdf = gdf.reset_index(drop=False)
            gdf['PLZ'] = gdf['index'].astype(str)
            gdf.drop(columns=['index'], inplace=True)

    return gdf


def load_germany_from_natural_earth() -> object:
    """Чтение 10m Natural Earth стран и возврат геометрии Германии (shapely)"""
    # Использование пути shapereader cartopy и geopandas для загрузки 10m admin_0_countries shapefile
    shp_path = shpreader.natural_earth(resolution='10m', category='cultural', name='admin_0_countries')
    countries = gpd.read_file(shp_path)
    # пробуем распространённые имена ISO-столбцов для сопоставления Германии
    iso_cols = [c for c in ('ISO_A3', 'ISO3', 'adm0_a3', 'ADM0_A3', 'iso_a3') if c in countries.columns]
    if iso_cols:
        germany = countries[countries[iso_cols[0]] == 'DEU']
    else:
        # запасной вариант: пробуем сопоставление по NAME или SOVEREIGNT
        name_cols = [c for c in ('NAME_EN', 'NAME', 'SOVEREIGNT', 'ADMIN') if c in countries.columns]
        if name_cols:
            germany = countries[countries[name_cols[0]].str.contains('Germany', case=False, na=False)]
        else:
            raise RuntimeError('Не удалось определить ISO или столбец имени в странах Natural Earth')
    if germany.empty:
        raise RuntimeError('Геометрия Германии не найдена в странах Natural Earth')
    return germany.iloc[0].geometry


def filter_plz_in_germany(plz_gdf, germany_geometry):
    """Фильтрация точек PLZ, находящихся в пределах Германии"""
    plz_gdf = plz_gdf.to_crs("EPSG:3857")  # Web Mercator для вычисления расстояний
    germany_geometry = gpd.GeoSeries([germany_geometry], crs="EPSG:4326").to_crs("EPSG:3857")[0]

    # Создание буфера вокруг Германии для включения PLZ около границы
    buffered_germany = germany_geometry.buffer(50000)  # буфер 50 км

    # Фильтрация PLZ в пределах буферизованной Германии
    plz_in_germany = plz_gdf[plz_gdf.intersects(buffered_germany)]
    return plz_in_germany


def calculate_coverage(selected_plz, germany_geometry, radius_km=30):
    """Вычисление процента Германии, покрытого выбранными PLZ"""
    radius_m = radius_km * 1000
    coverage = unary_union(selected_plz.geometry.buffer(radius_m))
    germany_geometry = gpd.GeoSeries([germany_geometry], crs="EPSG:4326").to_crs("EPSG:3857")[0]
    covered_area = germany_geometry.intersection(coverage).area
    total_area = germany_geometry.area
    return covered_area / total_area * 100


def _generate_grid_points_in_polygon(polygon_3857, step_m: float) -> np.ndarray:
    """Генерация регулярной сетки точек (x, y) внутри данного полигона (EPSG:3857)."""
    xmin, ymin, xmax, ymax = polygon_3857.bounds
    xs = np.arange(xmin, xmax + step_m, step_m)
    ys = np.arange(ymin, ymax + step_m, step_m)
    points = []
    # Простой цикл подходит для ~20k-60k точек
    for x in xs:
        for y in ys:
            p = Point(x, y)
            if polygon_3857.contains(p):
                points.append((x, y))
    if not points:
        return np.empty((0, 2), dtype=float)
    return np.array(points, dtype=float)


def select_plz_max_coverage(
    plz_gdf: gpd.GeoDataFrame,
    germany_geometry,
    radius_km: float = 30,
    coverage_target: float = 0.99,
    sample_size: int = 300,
    grid_step_km: float = 5.0,
    exclusion_factor: float = 0.75,
    max_steps: int = 10000,
    patience: int = 500,
    validate_every: int = 50,
    random_state: int = 42,
    verbose: bool = True
) -> gpd.GeoDataFrame:
    """
    Единый алгоритм выбора точек PLZ с максимизацией покрытой площади (приближённо).

    Использует регулярную сетку точек внутри Германии (EPSG:3857) как proxy площади.
    Жадный рандомизированный: на каждом шаге выбирается подмножество оставшихся кандидатов
    и выбирается тот, который покрывает больше всего непокрытых точек сетки в радиусе.

    Улучшения:
    - Периодические точные проверки покрытия (unary_union) для избежания преждевременной остановки.
    - Радиус исключения для подавления кластерных выборов: после выбора центра,
      отбрасываются все оставшиеся PLZ в пределах (exclusion_factor * radius_km).

    Возвращает: GeoDataFrame (EPSG:3857) выбранных точек PLZ.
    """
    rng = np.random.default_rng(random_state)

    # Гарантия проецированной CRS для расстояния/площади и подготовка геометрии
    work_gdf = plz_gdf.to_crs("EPSG:3857").copy()
    germany_3857 = gpd.GeoSeries([germany_geometry], crs="EPSG:4326").to_crs("EPSG:3857")[0]

    radius_m = float(radius_km) * 1000.0
    exclusion_m = float(exclusion_factor) * radius_m
    grid_step_m = float(grid_step_km) * 1000.0

    # 1) Построение точек сетки внутри Германии как proxy покрытия
    if verbose:
        print(f"Построение сетки {grid_step_km:.1f} км над Германией для аппроксимации покрытия...")
    grid_xy = _generate_grid_points_in_polygon(germany_3857, step_m=grid_step_m)
    if grid_xy.shape[0] == 0:
        raise RuntimeError("Не удалось сгенерировать точки сетки внутри полигона Германии.")
    if verbose:
        print(f"Точки сетки: {grid_xy.shape[0]:,}")

    # 2) KDTree по сетке для быстрых запросов в радиусе
    tree_grid = KDTree(grid_xy, leaf_size=40, metric='euclidean')

    # 3) Предварительное вычисление индексов соседних точек сетки для каждого кандидата PLZ
    plz_xy = np.column_stack([work_gdf.geometry.x.values, work_gdf.geometry.y.values])
    if verbose:
        print("Предварительное вычисление окрестностей покрытия кандидат-сетка...")
    neighbor_lists = tree_grid.query_radius(plz_xy, r=radius_m, return_distance=False)

    # 3b) KDTree по кандидатам PLZ для фильтрации исключений
    tree_plz = KDTree(plz_xy, leaf_size=40, metric='euclidean')

    # 4) Жадный рандомизированный цикл выбора по маске непокрытой сетки
    uncovered = np.ones(grid_xy.shape[0], dtype=bool)
    selected_mask = np.zeros(len(work_gdf), dtype=bool)
    active_mask = np.ones(len(work_gdf), dtype=bool)  # кандидаты всё ещё доступные
    no_improve_steps = 0

    pbar = tqdm(total=max_steps, disable=not verbose, desc="Выбор центров (макс-покрытие)")
    for step in range(max_steps):
        pbar.update(1)

        # Периодическая точная проверка покрытия
        if step % validate_every == 0 and selected_mask.any():
            selected_tmp = work_gdf.loc[selected_mask]
            exact_cov = calculate_coverage(selected_tmp, germany_geometry, radius_km) / 100.0
            if step % (validate_every * 2) == 0:
                pbar.set_postfix_str(f"exact~{exact_cov*100:.1f}%")
            if exact_cov >= coverage_target:
                break

        # Приближённая доля покрытия по сетке
        covered_frac = 1.0 - (uncovered.sum() / uncovered.size)
        if covered_frac >= coverage_target and selected_mask.any():
            selected_tmp = work_gdf.loc[selected_mask]
            exact_cov = calculate_coverage(selected_tmp, germany_geometry, radius_km) / 100.0
            if exact_cov >= coverage_target:
                break

        remaining_idx = np.flatnonzero(active_mask)
        if remaining_idx.size == 0:
            break

        # Выборка кандидатов из активного множества
        k = int(min(sample_size, remaining_idx.size))
        sample = rng.choice(remaining_idx, size=k, replace=False) if k > 0 else []

        # Оценка предельного выигрыша по сетке
        best_gain = 0
        best_idx = None
        for idx in sample:
            neigh = neighbor_lists[idx]
            if neigh.size == 0:
                continue
            gain = int(uncovered[neigh].sum())
            if gain > best_gain:
                best_gain = gain
                best_idx = idx

        if best_idx is None or best_gain == 0:
            no_improve_steps += 1
            if no_improve_steps >= patience:
                break
            else:
                continue

        # Принятие лучшего кандидата
        selected_mask[best_idx] = True
        # Пометка точек сетки как покрытых
        uncovered[neighbor_lists[best_idx]] = False
        # Исключение близлежащих PLZ для уменьшения кластеризации
        if exclusion_m > 0:
            near = tree_plz.query_radius(plz_xy[[best_idx]], r=exclusion_m, return_distance=False)[0]
            active_mask[near] = False
        # Также удаление выбранного индекса
        active_mask[best_idx] = False
        no_improve_steps = 0

    pbar.close()

    selected = work_gdf.loc[selected_mask].copy()
    selected = selected.set_crs("EPSG:3857")
    return selected


def visualize_coverage(selected_plz, germany_geometry, radius_km=30):
    """Визуализация покрытия с использованием Cartopy"""
    # Конвертация в WGS84 для визуализации
    selected_plz = selected_plz.to_crs("EPSG:4326")
    germany_geometry = gpd.GeoSeries([germany_geometry], crs="EPSG:4326").to_crs("EPSG:4326")[0]

    # Создание фигуры с проекцией Cartopy
    fig = plt.figure(figsize=(12, 12))
    ax = fig.add_subplot(1, 1, 1, projection=ccrs.EqualEarth())

    # Добавление элементов карты
    ax.add_feature(cfeature.LAND, facecolor='lightgray')
    ax.add_feature(cfeature.OCEAN, facecolor='lightblue')
    ax.add_feature(cfeature.COASTLINE, linewidth=0.5)
    ax.add_feature(cfeature.BORDERS, linestyle=':', linewidth=0.5)

    # Создание ShapelyFeature для Германии
    germany_feature = ShapelyFeature([germany_geometry], ccrs.PlateCarree(), facecolor='none', edgecolor='red', linewidth=1)
    ax.add_feature(germany_feature)

    # Создание и отрисовка кругов (теперь ВСЕ PLZ)
    radius_m = radius_km * 1000
    buffers_3857 = gpd.GeoSeries(selected_plz.geometry).to_crs("EPSG:3857").buffer(radius_m)
    buffers_4326 = gpd.GeoSeries(buffers_3857, crs="EPSG:3857").to_crs("EPSG:4326")

    # Отрисовка каждого буфера по отдельности, чтобы перекрытия отображались как слоистые прозрачности
    for geom in buffers_4326:
        ax.add_geometries([geom], crs=ccrs.PlateCarree(), facecolor='blue', edgecolor='none', alpha=0.2, zorder=2)

    # Отрисовка выбранных точек PLZ поверх
    ax.scatter(selected_plz.geometry.x, selected_plz.geometry.y, color='red', s=5, transform=ccrs.PlateCarree(), label='Выбранные PLZ', zorder=3)

    # Установка охвата для Германии с отступами
    ax.set_extent([4, 16, 47, 56], crs=ccrs.PlateCarree())

    # Добавление линий сетки
    gl = ax.gridlines(draw_labels=True, linestyle='--')
    gl.top_labels = False
    gl.right_labels = False

    plt.title(f"Покрытие PLZ Германии с {len(selected_plz)} центрами (радиус {radius_km} км)")
    plt.legend(loc='upper right')
    plt.show()


def run_analysis(plz_file_path, radius_km=30, coverage_target=0.99):
    """Главная функция для анализа с использованием единого алгоритма макс-покрытия"""
    # Загрузка геометрии Германии из Natural Earth
    print("Загрузка геометрии Германии из Natural Earth...")
    germany_geometry = load_germany_from_natural_earth()

    # Загрузка и обработка данных PLZ
    print("Загрузка данных PLZ...")
    plz_df = load_plz_data(plz_file_path)
    print(f"Загружено {len(plz_df)} немецких почтовых индексов")

    plz_gdf = create_plz_points(plz_df)

    print("Фильтрация PLZ в пределах Германии...")
    plz_in_germany = filter_plz_in_germany(plz_gdf, germany_geometry)
    print(f"Найдено {len(plz_in_germany)} PLZ в пределах или около Германии")

    # Выбор центров PLZ с использованием единого алгоритма
    print("Выбор центров PLZ с алгоритмом макс-покрытия...")
    selected_plz = select_plz_max_coverage(
        plz_in_germany, germany_geometry, radius_km=radius_km, coverage_target=coverage_target, sample_size=300, grid_step_km=2.0, exclusion_factor=0.75, max_steps=10000, patience=500, validate_every=50, random_state=42, verbose=True
    )

    # Вычисление покрытия (точное, единое объединение)
    coverage_percent = calculate_coverage(selected_plz, germany_geometry, radius_km)
    print(f"\nВыбрано {len(selected_plz)} центров PLZ, покрывающих {coverage_percent:.1f}% Германии")

    # Визуализация покрытия
    visualize_coverage(selected_plz, germany_geometry=germany_geometry, radius_km=radius_km)

    # Защитный возврат результатов: включаются только существующие столбцы
    sel = selected_plz.to_crs("EPSG:4326")
    desired = ['PLZ', 'lat', 'lon', 'geometry']
    available = [c for c in desired if c in sel.columns or c == 'geometry']
    # Гарантия наличия геометрии
    if 'geometry' not in available:
        available.append('geometry')

    return {
        'selected_plz': sel[available],
        'coverage_percent': coverage_percent,
        'total_plz_used': len(selected_plz)
    }

Этот код подходит для Jupyter Notebooks.

Пример использования

example-usage.py

results = run_analysis(
    plz_file_path="postleitzahlen.parquet",
    radius_km=35,
    coverage_target=0.995
)

results = run_analysis(
    plz_file_path="postleitzahlen.parquet",
    radius_km=35,
    coverage_target=0.995
)

Как найти ближайшие результирующие коды PLZ к заданному коду

Это опциональное дополнение принимает эталонный код PLZ некоторого места и генерирует DataFrame из следующих 50 ближайших выбранных PLZ. Это предполагает, что вы использовали Jupyter Notebook для запуска предыдущего кода.

Nearest PLZ.avif

find_nearest_plz.py

# Ячейка: поиск 50 выбранных PLZ, ближайших к заданному PLZ (пример: 63110)
plz_code = '10176'

# Гарантия наличия выбранных результатов; если нет, запускаем анализ один раз (лёгкая защита)
if 'results' not in globals():
    print('`results` не найден в пространстве имён ноутбука — запуск `run_analysis` (это может занять время)')
    results = run_analysis(plz_file_path="postleitzahlen.parquet", radius_km=35, coverage_target=0.995)

selected = results['selected_plz']

# Загрузка полной таблицы PLZ для нахождения координат эталонного PLZ
plz_df = load_plz_data("postleitzahlen.parquet")
plz_gdf = create_plz_points(plz_df)

# Точное совпадение по PLZ, иначе совпадение по префиксу
ref = plz_gdf[plz_gdf['PLZ'].astype(str) == str(plz_code)]
if ref.empty:
    ref = plz_gdf[plz_gdf['PLZ'].astype(str).str.startswith(str(plz_code))]

if ref.empty:
    raise ValueError(f"Эталонный PLZ {plz_code} не найден в наборе данных PLZ")

ref_point = ref.iloc[0].geometry

# Репроекция обоих в проецированную CRS для метрических расстояний (EPSG:3857)
selected_3857 = selected.to_crs("EPSG:3857").copy()
ref_3857 = gpd.GeoSeries([ref_point], crs="EPSG:4326").to_crs("EPSG:3857")[0]

# Вычисление расстояний (в метрах)
selected_3857['distance_m'] = selected_3857.geometry.distance(ref_3857)

# Если эталонный PLZ сам находится в выбранном множестве, исключаем его, чтобы получить ближайшие другие центры
if 'PLZ' in selected_3857.columns:
    selected_for_sort = selected_3857[selected_3857['PLZ'].astype(str) != str(plz_code)]
else:
    selected_for_sort = selected_3857

# Выбор 50 ближайших
nearest = selected_for_sort.nsmallest(50, 'distance_m').copy()

# Конвертация обратно в географические координаты и подготовка удобного DataFrame
nearest_geo = nearest.to_crs("EPSG:4326")
nearest_geo['lat'] = nearest_geo.geometry.y
nearest_geo['lon'] = nearest_geo.geometry.x
nearest_geo['distance_km'] = (nearest_geo['distance_m'] / 1000.0).round(3)

# Выбор столбцов для отображения: предпочтение человекочитаемым полям, если они есть
prefer_columns = ['PLZ', 'name', 'place', 'Ortsteil', 'stadtteil', 'place_name', 'lat', 'lon', 'distance_km']
available = [c for c in prefer_columns if c in nearest_geo.columns]
# если ни один из предпочтительных столбцов имени не существует, показываем все столбцы кроме geometry
if not available:
    available = [c for c in nearest_geo.columns if c != 'geometry']

nearest_50_df = nearest_geo[available].reset_index(drop=True)

print(f"Ближайшие {len(nearest_50_df)} выбранных PLZ к {plz_code} (топ {min(50, len(nearest_50_df))}):")
nearest_50_df.head(50)  # отображение DataFrame

# Ячейка: поиск 50 выбранных PLZ, ближайших к заданному PLZ (пример: 63110)
plz_code = '10176'

# Гарантия наличия выбранных результатов; если нет, запускаем анализ один раз (лёгкая защита)
if 'results' not in globals():
    print('`results` не найден в пространстве имён ноутбука — запуск `run_analysis` (это может занять время)')
    results = run_analysis(plz_file_path="postleitzahlen.parquet", radius_km=35, coverage_target=0.995)

selected = results['selected_plz']

# Загрузка полной таблицы PLZ для нахождения координат эталонного PLZ
plz_df = load_plz_data("postleitzahlen.parquet")
plz_gdf = create_plz_points(plz_df)

# Точное совпадение по PLZ, иначе совпадение по префиксу
ref = plz_gdf[plz_gdf['PLZ'].astype(str) == str(plz_code)]
if ref.empty:
    ref = plz_gdf[plz_gdf['PLZ'].astype(str).str.startswith(str(plz_code))]

if ref.empty:
    raise ValueError(f"Эталонный PLZ {plz_code} не найден в наборе данных PLZ")

ref_point = ref.iloc[0].geometry

# Репроекция обоих в проецированную CRS для метрических расстояний (EPSG:3857)
selected_3857 = selected.to_crs("EPSG:3857").copy()
ref_3857 = gpd.GeoSeries([ref_point], crs="EPSG:4326").to_crs("EPSG:3857")[0]

# Вычисление расстояний (в метрах)
selected_3857['distance_m'] = selected_3857.geometry.distance(ref_3857)

# Если эталонный PLZ сам находится в выбранном множестве, исключаем его, чтобы получить ближайшие другие центры
if 'PLZ' in selected_3857.columns:
    selected_for_sort = selected_3857[selected_3857['PLZ'].astype(str) != str(plz_code)]
else:
    selected_for_sort = selected_3857

# Выбор 50 ближайших
nearest = selected_for_sort.nsmallest(50, 'distance_m').copy()

# Конвертация обратно в географические координаты и подготовка удобного DataFrame
nearest_geo = nearest.to_crs("EPSG:4326")
nearest_geo['lat'] = nearest_geo.geometry.y
nearest_geo['lon'] = nearest_geo.geometry.x
nearest_geo['distance_km'] = (nearest_geo['distance_m'] / 1000.0).round(3)

# Выбор столбцов для отображения: предпочтение человекочитаемым полям, если они есть
prefer_columns = ['PLZ', 'name', 'place', 'Ortsteil', 'stadtteil', 'place_name', 'lat', 'lon', 'distance_km']
available = [c for c in prefer_columns if c in nearest_geo.columns]
# если ни один из предпочтительных столбцов имени не существует, показываем все столбцы кроме geometry
if not available:
    available = [c for c in nearest_geo.columns if c != 'geometry']

nearest_50_df = nearest_geo[available].reset_index(drop=True)

print(f"Ближайшие {len(nearest_50_df)} выбранных PLZ к {plz_code} (топ {min(50, len(nearest_50_df))}):")
nearest_50_df.head(50)  # отображение DataFrame

Check out similar posts by category: Algorithms, Geoinformatics

If this post helped you, please consider buying me a coffee or donating via PayPal to support research & publishing of new posts on TechOverflow

Buy me a coffee