Оценка равенства городской зелени с использованием открытого портала данных Вены

Assessing urban greenery equality using Vienna's open data portal.

Фото CHUTTERSNAP на Unsplash

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

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

Я буду исследовать два разных пространственных разрешения города — районы и переписные районы с использованием Esri Shapefiles, предоставленных Порталом открытых данных австрийского правительства. Я также буду включать табличные статистические данные (население и доход) в геопривязанные административные области. Затем я наложу административные области на официальный набор данных о зеленых зонах, записывая местоположение каждого зеленого пространства в геопространственном формате. Затем я объединю эту информацию и количественно оценю общую площадь зеленого пространства на душу населения каждого городского района. Наконец, я свяжу финансовое положение каждой области, зафиксированное годовым чистым доходом, с соотношением площади зеленого пространства на душу населения, чтобы определить, появляются ли какие-либо закономерности.

1. Источник данных

Давайте взглянем на Портал открытых данных австрийского правительства здесь.

Когда я писал эту статью, английский перевод веб-сайта не работал, поэтому вместо того, чтобы полагаться на мои давно забытые 12 летние уроки немецкого языка, я использовал DeepL для навигации по подстраницам и тысячам наборов данных.

Затем я собрал несколько файлов данных, как геопривязанных (Esri shapefiles), так и простых табличных данных, которые я буду использовать для последующего анализа. Собранные мной данные:

Границы — административные границы следующих пространственных единиц в Вене:

  • Административные границы Вены
  • Административные границы 23 районов Вены
  • Административные границы 250 переписных районов Вены

Землепользование — информация о местоположении зеленых зон и застроенных территорий:

  • Зеленый пояс Вены — город Вена, визуализация существующих и выделенных зеленых зон, состоящих из 1539 геопространственных полигональных файлов, ограничивающих зеленые пространства

Статистика — данные о населении и доходах, соответствующие социально-экономическому уровню территории:

  • Население по районам, ежегодно записываемое с 2002 года и сохраняемое разделенным на группы в возрасте 5 лет, пол и первоначальное гражданство
  • Население по переписным районам, ежегодно записываемое с 2008 года и сохраняемое разделенным на три нерегулярные возрастные группы, пол и происхождение
  • Средний чистый доход с 2002 года в районах Вены, выраженный в евро на сотрудника в год

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

2. Основное исследование данных

2.1 Административные границы

Сначала прочтите и визуализируйте различные файлы формата shape, содержащие каждый уровень административных границ, чтобы более подробно ознакомиться с городом:

папка = 'data'админ_город = gpd.read_file(папка + '/LANDESGRENZEOGD')админ_район = gpd.read_file(папка + '/BEZIRKSGRENZEOGD')админ_перепись = gpd.read_file(папка + '/ZAEHLBEZIRKOGD')display(админ_город.head(1))display(админ_район.head(1))display(админ_перепись.head(1))

Здесь мы замечаем, что имена столбцов BEZNR и ZBEZ соответствуют идентификаторам района и идентификаторам переписных районов соответственно. Неожиданно они хранятся/разбираются в разных форматах, numpy.float64 и str:

print(type(admin_district.BEZNR.iloc[0]))print(type(admin_census.ZBEZ.iloc[0]))pyth

Убедитесь, что у нас действительно есть 23 района и 250 переписных районов, как заявлено в документации к файлам данных:

print(len(set(admin_district.BEZNR)))print(len(set(admin_census.ZBEZ)))

Теперь визуализируйте границы — сначала города, затем его районов, а затем еще меньших переписных районов.

f, ax = plt.subplots(1,3,figsize=(15,5))admin_city.plot(ax=ax[0],       edgecolor = 'k',   linewidth = 0.5,   alpha = 0.9,   cmap = 'Reds')admin_district.plot(ax=ax[1],   edgecolor = 'k',   linewidth = 0.5,   alpha = 0.9,   cmap = 'Blues')admin_census.plot(ax=ax[2],     edgecolor = 'k',   linewidth = 0.5,   alpha = 0.9,   cmap = 'Purples')ax[0].set_title('Границы города')ax[1].set_title('Границы районов')ax[2].set_title('Границы переписных районов')

Этот код выводит следующую визуализацию Вены:

Различные административные уровни Вены. Изображение автора.

2.2 Зеленые зоны

Теперь давайте также посмотрим на распределение зеленых зон:

gdf_green  = gpd.read_file(folder + '/GRUENFREIFLOGD_GRUENGEWOGD')display(gdf_green.head(3))

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

Теперь визуализируйте это:

f, ax = plt.subplots(1,1,figsize=(7,5))gdf_green.plot(ax=ax,    edgecolor = 'k',   linewidth = 0.5,   alpha = 0.9,   cmap = 'Greens')ax.set_title('Зеленые зоны в Вене')

Этот код показывает, где находятся зеленые зоны в Вене:

Официальный зеленый пояс Вены. Изображение автора.

Можно отметить, что лесные сегменты все еще находятся в пределах административной границы, что означает, что не каждая часть города является застроенной и значительно населенной. Позже мы вернемся к этому при оценке зеленой площади на душу населения.

2.3 Статистические данные — население, доход

Наконец, давайте посмотрим на файлы статистических данных. Первое основное отличие состоит в том, что они не имеют географической привязки, а представляют собой простые таблицы csv:

df_pop_distr = pd.read_csv('vie-bez-pop-sex-age5-stk-ori-geo4-2002f.csv',   sep = ';',  encoding='unicode_escape',   skiprows = 1)df_pop_cens  = pd.read_csv('vie-zbz-pop-sex-agr3-stk-ori-geo2-2008f.csv',   sep = ';',  encoding='unicode_escape',   skiprows = 1)df_inc_distr = pd.read_csv('vie-bez-biz-ecn-inc-sex-2002f.csv',   sep = ';',  encoding='unicode_escape',   skiprows = 1)display(df_pop_distr.head(1))display(df_pop_cens.head(1))display(df_inc_distr.head(1))

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

3.1 Подготовка файлов статистических данных

Предыдущий подраздел показывает, что таблицы статистических данных используют разные соглашения об именовании — они имеют идентификаторы DISTRICT_CODE и SUB_DISTRICT_CODE вместо таких вещей, как BEZNR и ZBEZ. Однако, изучив документацию каждого набора данных, становится понятно, что легко перейти от одного к другому, для чего я представляю две короткие функции в следующей ячейке. Я буду одновременно обрабатывать данные на уровне районов и переписных районов.

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

# эти функции преобразуют идентификаторы района и района переписи, чтобы они были совместимы с теми, что находятся в файлах формы
def transform_district_id(x):   return int(str(x)[1:3])
def transform_census_district_id(x):   return int(str(x)[1:5])
# выберите последний год набора данных
df_pop_distr_2 = df_pop_distr[df_pop_distr.REF_YEAR \  ==max(df_pop_distr.REF_YEAR)]
df_pop_cens_2  = df_pop_cens[df_pop_cens.REF_YEAR \  ==max(df_pop_cens.REF_YEAR)]
df_inc_distr_2 = df_inc_distr[df_inc_distr.REF_YEAR \  ==max(df_inc_distr.REF_YEAR)]
# преобразование идентификаторов района
df_pop_distr_2['district_id'] = \  df_pop_distr_2.DISTRICT_CODE.apply(transform_district_id)
df_pop_cens_2['census_district_id'] = \  df_pop_cens_2.SUB_DISTRICT_CODE.apply(transform_census_district_id)
df_inc_distr_2['district_id'] = \  df_inc_distr_2.DISTRICT_CODE.apply(transform_district_id)
# агрегирование значений населения
df_pop_distr_2 = df_pop_distr_2.groupby(by = 'district_id').sum()
df_pop_distr_2['district_population'] = df_pop_distr_2.AUT + \  df_pop_distr_2.EEA  + df_pop_distr_2.REU  + df_pop_distr_2.TCN
df_pop_distr_2 = df_pop_distr_2[['district_population']]
df_pop_cens_2 = df_pop_cens_2.groupby(by = 'census_district_id').sum()
df_pop_cens_2['census_district_population'] = df_pop_cens_2.AUT \  + df_pop_cens_2.FOR
df_pop_cens_2 = df_pop_cens_2[['census_district_population']]
df_inc_distr_2['district_average_income'] = \  1000*df_inc_distr_2[['INC_TOT_VALUE']]
df_inc_distr_2 = \  df_inc_distr_2.set_index('district_id')[['district_average_income']]
# отображение окончательных таблиц
display(df_pop_distr_2.head(3))
display(df_pop_cens_2.head(3))
display(df_inc_distr_2.head(3))
# и унификация соглашений о названиях
admin_district['district_id'] = admin_district.BEZNR.astype(int)
admin_census['census_district_id'] = admin_census.ZBEZ.astype(int)
print(len(set(admin_census.ZBEZ)))

Проверьте вычисленные значения общей численности населения на двух уровнях агрегации:

print(sum(df_pop_distr_2.district_population))print(sum(df_pop_cens_2.census_district_population))

Оба этих значения должны дать одинаковый результат – 1931593 человека.

3.1. Подготовка файлов геопространственных данных

Теперь, когда мы закончили основную подготовку данных статистических файлов, пришло время сопоставить полигоны зеленой зоны с полигонами административной зоны. Затем давайте вычислим общую площадь зеленой зоны каждой административной зоны. Кроме того, я добавлю относительную площадь зеленой зоны каждой административной зоны из любопытства.

Чтобы получить значения в единицах СИ, нам нужно перейти к так называемой локальной системе координат, которая в случае Вены является EPSG:31282. Вы можете прочитать больше об этой теме, проекции карт и системы координат здесь и здесь.

# преобразование всех GeoDataFrames в локальную систему координат
admin_district_2 = \  admin_district[['district_id', 'geometry']].to_crs(31282)
admin_census_2 = \   admin_census[['census_district_id', 'geometry']].to_crs(31282)
gdf_green_2      = gdf_green.to_crs(31282)

Вычислите площадь административной единицы, измеренную в единицах СИ:

admin_district_2['admin_area'] = \  admin_district_2.geometry.apply(lambda g: g.area)
admin_census_2['admin_area'] = \   admin_census_2.geometry.apply(lambda g: g.area)
display(admin_district_2.head(1))
display(admin_census_2.head(1))

4. Вычисление отношения зеленой площади на душу населения

4.1 Вычисление покрытия зеленой площади в каждой административной единице

Для наложения двух геоданных GeoDataFrame административных границ и GeoDataFrame, содержащего полигоны зеленых зон, я буду использовать функцию overlay из GeoPandas. Затем я вычисляю площадь каждой зеленой зоны, попадающей в различные административные регионы. Далее я суммирую эти площади для каждой административной области, включая районы и переписные районы. В конечном шаге, на каждой единице разрешения, я добавляю ранее вычисленные официальные административные площади и вычисляю отношение общей площади к зеленой площади для каждого района и переписного района.

gdf_green_mapped_distr = gpd.overlay(gdf_green_2, admin_district_2)gdf_green_mapped_distr['green_area'] = \  gdf_green_mapped_distr.geometry.apply(lambda g: g.area) gdf_green_mapped_distr = \  gdf_green_mapped_distr.groupby(by = 'district_id').sum()[['green_area']]gdf_green_mapped_distr = \  gpd.GeoDataFrame(admin_district_2.merge(gdf_green_mapped_distr, left_on = 'district_id', right_index = True))gdf_green_mapped_distr['green_ratio'] = \  gdf_green_mapped_distr.green_area / gdf_green_mapped_distr.admin_areagdf_green_mapped_distr.head(3)

gdf_green_mapped_cens = gpd.overlay(gdf_green_2, admin_census_2)gdf_green_mapped_cens['green_area'] = \  gdf_green_mapped_cens.geometry.apply(lambda g: g.area)gdf_green_mapped_cens = \  gdf_green_mapped_cens.groupby(by = 'census_district_id').sum()[['green_area']]gdf_green_mapped_cens = \  gpd.GeoDataFrame(admin_census_2.merge(gdf_green_mapped_cens, left_on = 'census_district_id', right_index = True))gdf_green_mapped_cens['green_ratio'] = gdf_green_mapped_cens.green_area / gdf_green_mapped_cens.admin_areagdf_green_mapped_cens.head(3)

Наконец, визуализируем отношение зеленой площади по районам и переписным районам! Результаты кажутся вполне логичными, с высоким уровнем зелени на периферии и намного более низким уровнем в центральных районах. Кроме того, 250 переписных районов ясно показывают более детальную, детализированную картину характеристик разных районов, предлагая более глубокие и локализованные идеи для градостроителей. С другой стороны, информация на уровне районов, с десятикратно меньшим количеством пространственных единиц, показывает общие средние значения.

f, ax = plt.subplots(1,2,figsize=(17,5))gdf_green_mapped_distr.plot(ax = ax[0],   column = 'green_ratio',   edgecolor = 'k',   linewidth = 0.5,   alpha = 0.9,    legend = True,  cmap = 'Greens')gdf_green_mapped_cens.plot(ax = ax[1],   column = 'green_ratio',   edgecolor = 'k',   linewidth = 0.5,   alpha = 0.9,    legend = True,  cmap = 'Greens')

Этот блок кода выводит следующие карты:

Эти две карты показывают отношение зеленой площади в каждом районе / переписном районе в Вене. Изображение автора.

4.2 Добавление информации о населении и доходах для каждой административной единицы

В конечном шаге этого раздела давайте отобразим статистические данные в административных областях. Напоминание: у нас есть данные о населении как на уровне районов, так и на уровне переписных районов. Однако, я смог найти информацию о доходах (индикатор социально-экономического уровня) только на уровне районов. Это обычный компромисс в геопространственной науке. В то время как одно измерение (зелень) намного более информативно на более высоком разрешении (переписные районы), ограничения данных могут заставить нас все равно использовать более низкое разрешение.

display(admin_census_2.head(2))display(df_pop_cens_2.head(2))

gdf_pop_mapped_distr  = admin_district_2.merge(df_pop_distr_2, \  left_on = 'district_id', right_index = True)gdf_pop_mapped_cens  = admin_census_2.merge(df_pop_cens_2, \  left_on = 'census_district_id', right_index = True)gdf_inc_mapped_distr = admin_district_2.merge(df_inc_distr_2, \  left_on = 'district_id', right_index = True)f, ax = plt.subplots(1,3,figsize=(15,5))gdf_pop_mapped_distr.plot(column = 'district_population', ax=ax[0],     \  edgecolor = 'k', linewidth = 0.5, alpha = 0.9, cmap = 'Blues')gdf_pop_mapped_cens.plot(column = 'census_district_population', ax=ax[1], \  edgecolor = 'k', linewidth = 0.5, alpha = 0.9, cmap = 'Blues')gdf_inc_mapped_distr.plot(column = 'district_average_income', ax=ax[2],   \  edgecolor = 'k', linewidth = 0.5, alpha = 0.9, cmap = 'Purples')ax[0].set_title('district_population')ax[1].set_title('census_district_population')ax[2].set_title('district_average_incomee')

Этот блок кода приводит к следующей фигуре:

Различная статистическая информация о районах Вены. Изображение автора.

4.3. Вычисление площади зеленой зоны на душу населения

Давайте подведем итоги того, что у нас есть сейчас, все объединено в достойные shape-файлы, соответствующие районам и переписным районам Вены:

На уровне районов у нас есть соотношение зеленой зоны, население и данные о доходе

На уровне переписных районов у нас есть соотношение зеленой зоны и население

Чтобы просто учесть зеленое равенство, я объединяю информацию о абсолютном размере зеленой зоны и населении в районах и переписных районах и вычисляю общую площадь зеленой зоны на душу населения.

Давайте посмотрим на наш ввод – охват зеленой зоны и население:

# график для районовf, ax = plt.subplots(1,2,figsize=(10,5))gdf_green_mapped_distr.plot(  ax = ax[0],   column = 'green_ratio',   edgecolor = 'k',   linewidth = 0.5,   alpha = 0.9,    cmap = 'Greens')gdf_pop_mapped_distr.plot(  ax = ax[1],   column = 'district_population',    edgecolor = 'k',   linewidth = 0.5,   alpha = 0.9,    cmap = 'Reds')ax[0].set_title('green_ratio')ax[1].set_title('district_population')# график для переписных районовf, ax = plt.subplots(1,2,figsize=(10,5))gdf_green_mapped_cens.plot(  ax = ax[0],   column = 'green_ratio',   edgecolor = 'k',   linewidth = 0.5,   alpha = 0.9,    cmap = 'Greens')gdf_pop_mapped_cens.plot(ax = ax[1],   column = 'census_district_population',  edgecolor = 'k',   linewidth = 0.5,   alpha = 0.9,    cmap = 'Reds')ax[0].set_title('green_ratio')ax[1].set_title('district_population')

Этот блок кода приводит к следующей фигуре:

Уровень зеленой зоны и население в Вене на уровне районов и переписных районов. Изображение автора.

Для вычисления площади зеленой зоны на душу населения я сначала объединю фреймы данных о зеленых насаждениях и населении в следующих шагах. Я сделаю это на примере переписных районов, потому что их более высокое пространственное разрешение позволяет нам лучше наблюдать возникающие закономерности (если таковые имеются). Убедитесь, что мы не делим на ноль и также следуем здравому смыслу; давайте удалим те области, которые не имеют населения.

gdf_green_pop_cens = \  gdf_green_mapped_cens.merge(gdf_pop_mapped_cens.drop( \    columns = ['geometry', 'admin_area']), left_on = 'census_district_id',\    right_on = 'census_district_id')[['census_district_id', \    'green_area', 'census_district_population',  'geometry']]gdf_green_pop_cens['green_area_per_capita'] = \  gdf_green_pop_cens['green_area'] / \  gdf_green_pop_cens['census_district_population']gdf_green_pop_cens = \  gdf_green_pop_cens[gdf_green_pop_cens['census_district_population']>0]f, ax = plt.subplots(1,1,figsize=(10,7))gdf_green_pop_cens.plot(  column = 'green_area_per_capita',   ax=ax,   cmap = 'RdYlGn',   edgecolor = 'k',   linewidth = 0.5)admin_district.to_crs(31282).plot(\  ax=ax, color = 'none', edgecolor = 'k', linewidth = 2.5)

Этот блок кода приводит к следующей фигуре:

Каждый район перекрашен в зависимости от оценки зеленой площади на душу населения. Изображение автора.

Давайте немного изменить визуализацию:

f, ax = plt.subplots(1,1,figsize=(11,7))ax.set_title('Зеленая площадь на душу населения в районах переписи населения Вены',   fontsize = 18, pad = 30)gdf_green_pop_cens.plot(  column = 'green_area_per_capita',   ax=ax,   cmap = 'RdYlGn',   edgecolor = 'k',   linewidth = 0.5,   legend=True,   norm=matplotlib.colors.LogNorm(\    vmin=gdf_green_pop_cens.green_area_per_capita.min(), \    vmax=gdf_green_pop_cens.green_area_per_capita.max()), )admin_district.to_crs(31282).plot(  ax=ax, color = 'none', edgecolor = 'k', linewidth = 2.5)

Этот блок кода приводит к следующей диаграмме:

Каждый район перекрашен в зависимости от оценки зеленой площади на душу населения. Изображение автора.

То же самое для районов:

# вычисление оценок зеленой площади на душу населенияgdf_green_pop_distr = \  gdf_green_mapped_distr.merge(gdf_pop_mapped_distr.drop(columns = \  ['geometry', 'admin_area']), left_on = 'district_id', right_on = \   'district_id')[['district_id', 'green_area', 'district_population', \   'geometry']]gdf_green_popdistr = \  gdf_green_pop_distr[gdf_green_pop_distr.district_population>0]gdf_green_pop_distr['green_area_per_capita'] = \  gdf_green_pop_distr['green_area'] / \  gdf_green_pop_distr['district_population']# визуализация карты на уровне районовf, ax = plt.subplots(1,1,figsize=(10,8))ax.set_title('Зеленая площадь на душу населения в районах Вены', \   fontsize = 18, pad = 26)gdf_green_pop_distr.plot(column = 'green_area_per_capita', ax=ax, \  cmap = 'RdYlGn', edgecolor = 'k', linewidth = 0.5, legend=True, \  norm=matplotlib.colors.LogNorm(vmin=\  gdf_green_pop_cens.green_area_per_capita.min(), \  vmax=gdf_green_pop_cens.green_area_per_capita.max()), )admin_city.to_crs(31282).plot(ax=ax, \  color = 'none', edgecolor = 'k', linewidth = 2.5)

Этот блок кода приводит к следующей диаграмме:

Каждый район перекрашен в зависимости от оценки зеленой площади на душу населения. Изображение автора.

Хотя значимые тенденции явны – наружный край, больше зеленого пространства для всех, встроенный центр, обратный. Тем не менее, эти две диаграммы, особенно более подробная на уровне районов переписи населения, ясно показывают различие в количестве зеленого пространства, которым люди наслаждаются в разных районах. Дальнейшие исследования и включение дополнительных источников данных, например, о использовании земли, могут помочь лучше объяснить, почему эти районы более зеленые или населенные. На данный момент давайте насладимся этой картой и надеемся, что каждый найдет нужное количество зелени дома!

# объединение данных о зеленых площадях, населении и финансовых данныхgdf_district_green_pip_inc = \  gdf_green_pop_distr.merge(gdf_inc_mapped_distr.drop(columns = \  ['geometry']))

Визуализация связи между финансовыми и зелеными измерениями:

f, ax = plt.subplots(1,1,figsize=(6,4))
ax.plot(gdf_district_green_pip_inc.district_average_income, \
         gdf_district_green_pip_inc.green_area_per_capita, 'o')
ax.set_xscale('log')
ax.set_yscale('log')
ax.set_xlabel('Средний доход района')
ax.set_ylabel('Зеленая площадь на душу населения')

Результат выполнения этого кода – следующая диаграмма рассеяния:

Построение отношения среднего чистого дохода и площади зеленых на душу населения в районах Вены. Изображение автора.

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

print(spearmanr(gdf_district_green_pip_inc.district_average_income, gdf_district_green_pip_inc.green_area_per_capita))
print(pearsonr(gdf_district_green_pip_inc.district_average_income, gdf_district_green_pip_inc.green_area_per_capita))

Из-за тяжелого распределения финансовых данных я бы отнесся более серьезно к коэффициенту Спирмена (0.13), но даже коэффициент Пирсона (0.30) указывает на относительно слабую тенденцию, что соответствует моим предыдущим наблюдениям.