В этом посте я покажу, как агрегировать информацию из двух разных GeoDataFrames (слоев) с помощью GeoPandas. Я покажу один из многих способов сделать это на примере кода. Но сначала давайте посмотрим на ситуацию, когда нам может понадобиться выполнить такую ​​операцию.

Obs 01: Блокнот и данные, использованные в этом посте, доступны на моей странице Github. Если хотите, можете подписаться на этот пост вместе с блокнотом.

Obs 02: Если вы хотите отобразить записную книжку в своем браузере, используйте эту ссылку.

Проблема

Допустим, вы работаете в Carrefour в городе Сан-Паулу и хотите знать, сколько людей живет в радиусе 1 км от всех магазинов города Сан-Паулу и общий ежемесячный доход. Имея эту информацию, вы хотите расставить приоритеты в маркетинговой кампании для магазинов с самым высоким ежемесячным доходом на душу населения.

У нас есть кадр данных, содержащий информацию о местонахождении всех магазинов в городе Сан-Паулу. Мы собираемся извлечь данные о численности населения и общем ежемесячном доходе из переписи населения Бразилии, подготовленной БИГС.

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

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

Как я упоминал выше, я собрал информацию о расположении магазинов Carrefour в городе Сан-Паулу.

# Importing libraries
import pandas as pd
import numpy as np
import geopandas as gpd
import folium

# Loading Carrefour Store data
df_carrefour = pd.read_excel(r'https://github.com/AlvaroMatsuda/medium_posts/blob/main/spatial%20aggregation/data/carrefour_stores_latlong.xlsx?raw=true')
df_carrefour.head()

Этот кадр данных представляет собой объект Pandas DataFrame, который мы преобразуем в объект GeoPandas GeoDataFrame. Для этого нам нужно извлечь WKT (Well Known Text), в котором хранится пространственная информация из столбцов lat и long. После этого мы создадим буфер радиусом 1 км для каждого магазина.

# Converting Pandas DataFrame into GeoPandas GeoDataFrame

# Construct geometry Column from lat and long coordinates using geopandas
geometry = gpd.points_from_xy(df_carrefour['long'], df_carrefour['lat'], crs="EPSG:4979")

# Adding new column with the geometry
df_carrefour['geometry'] = geometry

# Constructing GeoDataFrame
gdf_carrefour = gpd.GeoDataFrame(df_carrefour, geometry=df_carrefour['geometry'])
gdf_carrefour.head()

С объектом GeoDataFrame мы можем использовать метод explore() для построения карты с расположением магазинов Carrefour:

# Plotting the map to visualize
gdf_carrefour.explore()

Теперь нам нужно создать буфер в 1 км для каждого магазина.

# Creating radius buffer
# Converting CRS to a meter based CRS
gdf_carrefour.to_crs(crs=3857, inplace=True) 

# Creating 1km buffer column with WKT geometry
gdf_carrefour['buffer_geom'] = gdf_carrefour.buffer(1000.0) 

# Converting back to original CRS
gdf_carrefour.to_crs(crs=4979, inplace=True) 

# Setting the geometry column to the buffer geometry
gdf_carrefour.set_geometry("buffer_geom", inplace=True)

gdf_carrefour.head()

# Plotting the map to visualize
gdf_carrefour.explore()

Вот и вся подготовка к фрейму данных магазина Carrefour. Теперь нам нужно подготовить данные переписных участков.

Прежде чем мы начнем подготовку к данным переписи, я должен упомянуть, что данные, которые нам нужны для этой задачи, находятся в 3 разных файлах:

  1. Формы переписных участков (многоугольники)
  2. csv с информацией о подсчете населения
  3. csv файл с информацией о доходах

Для начала мы загружаем формы/многоугольники переписных участков:

# Loading Brazilian census tract data
census_tract_shapes = gpd.GeoDataFrame.from_file(r'/vsicurl/https://raw.githubusercontent.com/AlvaroMatsuda/medium_posts/main/spatial%20aggregation/data/setor_censitario_sp_2010.shp')

# Setting the CRS
census_tract_shapes = census_tract_shapes.to_crs('EPSG:4979')

# Filtering to get only shapes from São Paulo city
sp_census_tract_shapes = census_tract_shapes.loc[census_tract_shapes['NM_MUNICIP'] == 'SÃO PAULO']

sp_census_tract_shapes.head()

# Plottin the map to visualize
sp_census_tract_shapes.explore()

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

# Loading total Monthly Income dataframe
# The column that holds the total monthly income is: V003
total_monthly_income_df = pd.read_csv(r'https://github.com/AlvaroMatsuda/medium_posts/blob/main/spatial%20aggregation/data/DomicilioRenda_SP1.csv?raw=true', sep=';', encoding='latin1')

# Converting Column "Cod_setor" from int to str
total_monthly_income_df['Cod_setor'] = total_monthly_income_df['Cod_setor'].astype('str')

total_monthly_income_df.head()

# Loading population dataframe
# The column that holds the population count is: V002
population_df = pd.read_csv(r'https://github.com/AlvaroMatsuda/medium_posts/blob/main/spatial%20aggregation/data/Pessoa11_SP1.csv?raw=true', sep=';', encoding='latin1')

# Converting Column "Cod_setor" from int to str
population_df['Cod_setor'] = population_df['Cod_setor'].astype('str')

population_df.head()

Все фреймы данных имеют один общий столбец, который представляет собой код тракта «Cod_setor» и «CD_GEOCODI». Мы собираемся присоединиться к этим фреймам данных в этом столбце в качестве ключа для получения необходимой информации о фрейме данных, содержащем фигуры/многоугольники.

# Joining shape dataframe with population and total_monthly_income dataframes
sp_census_tract = (
    sp_census_tract_shapes
    .merge(population_df[['Cod_setor', 'V002']], left_on='CD_GEOCODI', right_on='Cod_setor', how='left') # Joining dataframe
    .merge(total_monthly_income_df[['Cod_setor', 'V003']], left_on='CD_GEOCODI', right_on='Cod_setor', how='left') # Joining dataframe
    .drop(columns=['Cod_setor_x', 'Cod_setor_y']) # Dropping duplicated columns
    .rename(columns={'V002': 'pop_count', 'V003': 'total_monthly_income'}) # Renaming columns
    )

# Renaming columns to be all in lowercase
sp_census_tract.columns = map(str.lower, sp_census_tract.columns)

# Convertin "X" values to 0
sp_census_tract['pop_count'] = np.where(sp_census_tract['pop_count'] == "X", 0, sp_census_tract['pop_count'])
sp_census_tract['total_monthly_income'] = np.where(sp_census_tract['total_monthly_income'] == "X", 0, sp_census_tract['total_monthly_income'])

# Converting STR to FLOAT of pop_count and total_monthly_income columns
sp_census_tract['pop_count'] = sp_census_tract['pop_count'].astype('float64')
sp_census_tract['total_monthly_income'] = sp_census_tract['total_monthly_income'].astype('float64')

sp_census_tract.head()

Мы добавили столбцы pop_count и total_montly_income во фрейм данных формы, который содержит столбец геометрии. Вот и все необходимые манипуляции.

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

Объединение информации двух GeoDataFrames

Прежде чем мы сделаем агрегацию для всех буферов, я сделаю агрегацию для одного буфера, чтобы было проще понять сделанные манипуляции и код.

Прежде всего, давайте получим один буфер:

# Selecting a single buffer
single_buffer = gdf_carrefour[gdf_carrefour['store_name'] == 'Market Paraíso']

# Plottin the map to visualize
single_buffer.explore()

Теперь получим все переписные участки, пересекающие буфер:

геометрические объекты пересекаются, если они имеют какую-либо общую границу или внутреннюю точку. (определение из стройной документации)

# Filtering census tract that intersects the buffer
intersected_census_tract = sp_census_tract[sp_census_tract['geometry'].intersects(single_buffer.loc[8, '1km_buffer_geom'])]

# Plottin the map to visualize
intersected_census_tract.explore()

Давайте визуализируем оба слоя вместе на одной карте:

# Creating a folium tile base map
m = folium.Map(location=[-23.568865965761308, -46.629872859184125], zoom_start=12, tiles='OpenStreetMap', crs='EPSG3857', width='60%', height='60%')

# Adding intesected census_tract layer to the map
intersected_census_tract.explore(m=m)

# Adding buffer layer to the map
single_buffer.explore(m=m,style_kwds={'fillOpacity': 0.0, 'color': '#CE2020'})

# Showing the map
m

На изображении выше показано, что мы выбрали все участки переписи (синие многоугольники), которые пересекают буфер (красный многоугольник). Теперь нам нужно просуммировать столбцы pop_count и total_monthly_income, чтобы получить необходимую информацию. Однако, прежде чем мы суммируем эти столбцы, мы должны подумать, что делать с участками, площадь которых находится за пределами буфера. Должны ли мы рассматривать только тракты, которые на 100% внутри? Должны ли мы рассматривать все эти трактаты? Или мы должны рассматривать только часть информации, хранящейся в трактах, имеющих площадь за пределами буфера?

Есть много способов, которыми мы можем работать с этими полигонами с областью вне буфера.

  • Мы могли бы использовать центроид полигона: если центроид находится внутри буфера, мы рассматриваем всю информацию, хранящуюся в этом полигоне.
  • Учитывайте только те полигоны, которые на 100 % находятся внутри буфера.
  • Рассмотрим часть данных, хранящихся в полигоне, в соответствии с долей площади внутри буфера.

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

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

# Calculating the proportion of tract area inside the buffer
intersected_census_tract['prop_area'] = intersected_census_tract['geometry'].intersection(single_buffer.loc[8, '1km_buffer_geom']).area / intersected_census_tract['geometry'].area

# Rounding number to be 1 (100%) 
intersected_census_tract['prop_area'] = np.where(intersected_census_tract['prop_area'] >= 0.998, 1 , intersected_census_tract['prop_area'])

intersected_census_tract.head()

Чтобы объяснить, что делает приведенный выше код, в основном он вычисляет площадь внутри буфера (как карта справа) и делит результат на общую площадь этого участка. После этого мы создаем новый столбец для хранения доли площади. Кроме того, код округляет число до 1, если значение в столбце prop_area больше 0,998. Мы делаем это, потому что, когда мы делаем деление, вычисленная площадь не является точной, имея некоторые незначительные различия, поэтому, даже если тракт находится внутри на 100%, значение будет 0,9999 из-за этой незначительной разницы.

После этого мы можем использовать столбец prop_area для расчета доли информации, которая будет учитываться при суммировании pop_count. и total_monthly_income.

# Calculating the proportion of the pop_count according to the prop_area
intersected_census_tract['prop_pop_count'] = round(intersected_census_tract['pop_count'] * intersected_census_tract['prop_area'])

# Calculating the proportion of the total_monthly_income according to the prop_area
intersected_census_tract['prop_total_monthly_income'] = round(intersected_census_tract['total_monthly_income'] * intersected_census_tract['prop_area'])

intersected_census_tract.head()

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

# Calculating the sum of pop_count and total_monthly_income inside the buffer
single_buffer['1km_buffer_pop_count'] = intersected_census_tract['prop_pop_count'].sum()
single_buffer['1km_buffer_total_monthly_income'] = intersected_census_tract['prop_total_monthly_income'].sum()

single_buffer

Теперь, когда мы сделали всю обработку для одного буфера, я покажу код для всех буферов.

# Iterating over each buffer
for index, value in gdf_carrefour.iterrows():

    # Selecting Tracts that intersect the buffer
    df_intersected = sp_census_tract[sp_census_tract['geometry'].intersects(value['1km_buffer_geom'])]

    # Calculating the proportion of the area inside the buffer
    df_intersected['prop_area'] = df_intersected['geometry'].intersection(value['1km_buffer_geom']).area / df_intersected['geometry'].area

    # Rounding number to be 1 (100%) 
    df_intersected['prop_area'] = np.where(df_intersected['prop_area'] >= 0.998, 1, df_intersected['prop_area'])

    # Calculating the proportion of the pop_count and total_monthly_income according to the prop_area
    df_intersected[['prop_pop_count', 'prop_total_monthly_income']] = df_intersected.loc[:, ['pop_count', 'total_monthly_income']].multiply(df_intersected.loc[:,'prop_area'], axis=0)

    # Calculating total pop_count and total monthly income inside the buffer
    gdf_carrefour.loc[index, ['1km_buffer_pop_count', '1km_buffer_total_monthly_income']] = list(df_intersected[['prop_pop_count', 'prop_total_monthly_income']].sum())
    
gdf_carrefour.head()

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

# Calculating monthly income per capita
gdf_carrefour['monthly_income_per_capita'] = gdf_carrefour['1km_buffer_total_monthly_income'] / gdf_carrefour['1km_buffer_pop_count']

# Sorting DataFrame to show 5 stores with highest monthly income per capita
gdf_carrefour.sort_values(by='monthly_income_per_capita', ascending=False).head()

Благодаря этому мы могли бы расставить приоритеты в маркетинговой кампании для этих 5 магазинов, чтобы привлечь новых клиентов и заставить клиентов тратить больше на эти магазины.

# Creating a copy to use as store point layer
point_gdf_carrefour = gdf_carrefour.copy()

# Setting the geometry to the "geometry" columns (geometry representing points of the location of stores)
point_gdf_carrefour = point_gdf_carrefour.set_geometry('geometry')

# Creating a folium tile base map
m = folium.Map(location=[-23.568865965761308, -46.629872859184125], zoom_start=12, tiles='OpenStreetMap', crs='EPSG3857', width='70%', height='70%')

# Adding buffer layer to the map
gdf_carrefour.sort_values(by='monthly_income_per_capita', ascending=False).head().explore(style_kwds={'fillOpacity': 0.0}, m=m)

# Adding Store location layer to the map
point_gdf_carrefour.loc[gdf_carrefour.sort_values(by='monthly_income_per_capita', ascending=False).head().index, :].explore(style_kwds={'color': '#000000', 'weight': 8}, m=m)

# Showing the map
m

Заключение

На рынке ГИС (географических информационных систем) такого рода манипуляции используются очень часто, и есть много других вариантов их использования.

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

Я надеюсь, что этот пост вам чем-то поможет, и не забывайте, что блокнот с кодом и данными, использованными в этом посте, находится на моем Github.



medium_posts/spatial aggregation at main · AlvaroMatsuda/medium_posts
Вы не можете выполнить это действие в данный момент. Вы вошли в другую вкладку или окно. Вы вышли на другой вкладке или…github.com



Обо мне

Я географ, в настоящее время работаю специалистом по данным. По этой причине я сторонник науки о пространственных данных.

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