1. Вариации на тему НРМ - Advanced Spatial Analysis

Неравномерно-районированная модель города была предложена Александром Аркадьевичем Высоковским и описывает развитие городской территории как изначально неоднородный процесс. Согласно этой концепции, рост и концентрация функций обычно тяготеют к основному центру города — чаще всего историческому ядру. В то же время на периферии возникают вторичные очаги городской активности: локальные центры, ориентированные на повседневные нужды жителей — торговлю, услуги, рабочие места и социальную инфраструктуру.

В этом разделе мы не будем подробно рассматривать модель. А лишь выполним одно небольшое “упражнение”, которое покажет один из возможных подходов выделения локальных центров на основе выборки данных.

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

Импортируем библиотеки

import geopandas as gpd

import osmnx as ox


from libpysal.weights import Queen
from pysal.lib import weights

Загружаем данные о плотности населения Краснодара, агрегированные по квадратной сетке регулярных ячеек

sq_grid = gpd.read_file('./data/krs_pop_sqgrid.geojson')

sq_grid.explore(tiles='cartodbpositron')

Записывем СК сетки в переменную и название города, с которым мы работаем

krs_crs = sq_grid.estimate_utm_crs()
place="Krasnodar city"

Загружаем точки торговли из osm

tags_retail = {
    'shop': [
        'supermarket',
        'convenience',
        'grocery'
    ]
}

retail = ox.features_from_place(place, tags_retail)

retail.explore(tiles='cartodbpositron')

Фильтруем данные и считаем количество точек по ячейкам регулярной сетки

retail_pts = retail[retail.geometry.type == "Point"].copy()

# убираем пустые геометрии
retail_pts = retail_pts[~retail_pts.geometry.is_empty & retail_pts.geometry.notna()].copy()

# перепроецируем в UTM-зону
retail_pts = retail_pts.to_crs(krs_crs)

# объединяем с ячейками сетки
join = gpd.sjoin(
    retail_pts[["geometry"]],
    sq_grid[["geometry"]],
    how="left",
    predicate="within"
)

# считаем количество точек по каждой ячейке
counts = join.groupby("index_right").size().rename("retail_cnt")

# записываем посчитанное в исходную сетку
sq_grid["retail_cnt"] = sq_grid.index.map(counts)

# смотрим на результат
sq_grid.explore(column="retail_cnt", tiles='cartodbpositron',
        missing_kwds={
        "color": "lightgrey",
        "label": "No data"
    })

Оставляем только ячейки внутри Краснодара

krs_border = ox.geocode_to_gdf(place)
krs_border_utm = krs_border.to_crs(krs_crs)

krs_grid = sq_grid[sq_grid.intersects(krs_border_utm.geometry.iloc[0])].copy()

krs_grid.explore(column="retail_cnt", tiles='cartodbpositron')

1. Выделение локальных центров¶

1.1 Расчёт пространственного лага (среднего значения в соседних ячейках)¶

df = krs_grid.copy()

#изучаемое значение
x_col = "retail_cnt"  

# убираем Nan
df[x_col] = df[x_col].fillna(0)

# веса Queen
w = Queen.from_dataframe(df)

print("Number of islands (no neighbors):", len(w.islands))

#row-стандартизация весов
w.transform = "r"

# пространственный лаг 
x = df[x_col].values
wx = weights.lag_spatial(w, x)# Wx_i

df["W_" + x_col] = wx

Number of islands (no neighbors): 0

/var/folders/ry/9bb7wrz54vq_kn2ytlj6ynzm0000gn/T/ipykernel_63116/907929596.py:10: FutureWarning: `use_index` defaults to False but will default to True in future. Set True/False directly to control this behavior and silence this warning
  w = Queen.from_dataframe(df)

1.2 Вычсление отклонений от пространственного лага¶

df["dev"] = df[x_col] - df["W_" + x_col]    # d_i = x_i - Wx_i

df.explore(
    column="dev",
    legend=True,
    tiles="cartodbpositron"
)

1.3 Определение центров¶

# стандартизация (z-score)
dev_mean = df["dev"].mean()
dev_std = df["dev"].std(ddof=0)
df["dev_z"] = (df["dev"] - dev_mean) / dev_std

# центры по порогу
df["center_1sd"] = df["dev_z"] >= 1
df["center_2sd"] = df["dev_z"] >= 2

df[["dev", "dev_z", "center_1sd", "center_2sd"]].head()

1.4 Изучаем результат¶

df["center_type"] = "no center"
df.loc[df["dev_z"] >= 1, "center_type"] = "center (≥1σ)"
df.loc[df["dev_z"] >= 2, "center_type"] = "core center (≥2σ)"

df.explore(
    column="center_type",
    categorical=True,
    legend=True,
    tiles="cartodbpositron",
    tooltip=[x_col, "dev_z"]
)