<jupyter_start><jupyter_text>Study area[Natal, Brazil](https://en.wikipedia.org/wiki/Natal,_Rio_Grande_do_Norte)<jupyter_code>city = api.get('(-8.02, -41.01, -3.0, -33.0)',
                 largest_geom=True,
                 name='Natal')
city.plot();<jupyter_output>Requesting name=Natal
Geohunter: [TIMELOG] get -- {'largest_geom': True, 'name': 'Natal'} -- Completed in 10.0411s<jupyter_text>Grids<jupyter_code>resolution = 1
p_grid = gh.util.make_gridpoints(city, resolution=resolution/2)
p_grid = p_grid.drop(columns=['lat', 'lon'])
p_grid.plot(markersize=1);
s_grid = gh.util.make_gridsquares(city, resolution)[['geometry']]
s_grid.plot(edgecolor='black')<jupyter_output><empty_output><jupyter_text>Extracting Points-of-Interest data[OpenStreetMap's map features](https://wiki.openstreetmap.org/wiki/Map_Features)<jupyter_code>poi = api.get(city, 
    amenity=['school', 'hospital', 'restaurant', 'police', 'place_of_worship'],
    leisure='*',
    highway=['primary', 'residential', 'bus_stop'],
    tourism='*',      
    natural=['sand', 'wood', 'beach'],
    shop='*'
)
poi.loc[poi['key'].isin(['leisure', 'tourism', 'shop']), 'item'] = '*'
poi.plot(column='key')
poi.head()
poi.groupby(['key','item']).size()<jupyter_output><empty_output><jupyter_text>Extracting geographic features/covariates  1. Quadrat method<jupyter_code>x = []
for pair in poi[['key', 'item']].drop_duplicates().iterrows():
    key, item = pair[1].values
    x.append(gpd.sjoin(poi.loc[(poi['key']==key) & (poi['item']==item)], s_grid)\
            .groupby('index_right').size().to_frame(f'{key}_{item}'))
    
X_quadrat = pd.concat(x, axis=1)\
                    .join(s_grid, how='outer')\
                    .fillna(0)
X_quadrat = gpd.GeoDataFrame(X_quadrat)
feature = 'amenity_school'
X_quadrat.plot(column=feature)<jupyter_output><empty_output><jupyter_text>2. KDE method<jupyter_code>area_sgrid = gpd.GeoDataFrame(geometry=[s_grid.geometry.unary_union])
x = []
for pair in poi[['key', 'item']].drop_duplicates().iterrows():
    key, item = pair[1].values
    x.append(gh.util.kde_interpolation(poi.loc[(poi['key']==key) & (poi['item']==item)], 
                                       bw='silverman', area=area_sgrid, resolution=resolution/2)['density']\
                                       .to_frame(f'{key}_{item}'))
X_kde = pd.concat(x, axis=1)
X_kde = p_grid.join(X_kde)
feature = 'amenity_school'
X_kde.plot(column=feature)<jupyter_output><empty_output><jupyter_text>Average values within the square grid cells:<jupyter_code>X_kde = gpd.sjoin(X_kde, s_grid)\
                .groupby('index_right').mean()\
                .join(s_grid, how='outer')\
                .fillna(0)
X_kde = gpd.GeoDataFrame(X_kde)
X_kde.plot(feature)
X_kde.head()<jupyter_output><empty_output><jupyter_text>Qualitative assessment<jupyter_code>fig, ax = plt.subplots(ncols=3, nrows=3, figsize=(14,12))

for i in range(3):
    for j in range(3):
        s_grid.plot(edgecolor='black', color='white', ax=ax[i][j])
        

poi.loc[(poi['key']=='amenity') & (poi['item']=='hospital')]['geometry'].centroid.plot(ax=ax[0][0], color='red', markersize=5)
poi.loc[(poi['key']=='highway') & (poi['item']=='bus_stop')]['geometry'].centroid.plot(ax=ax[0][1], color='green', markersize=5)
poi.loc[(poi['key']=='highway') & (poi['item']=='residential')]['geometry'].centroid.plot(ax=ax[0][2], markersize=3)

X_quadrat.plot(column='amenity_hospital', ax=ax[1][0], legend=True, edgecolor='black', cmap='Reds')
X_quadrat.plot(column='highway_bus_stop', ax=ax[1][1], legend=True, edgecolor='black', cmap='Greens')
X_quadrat.plot(column='highway_residential', ax=ax[1][2], legend=True, edgecolor='black', cmap='Blues')

X_kde.plot(column='amenity_hospital', ax=ax[2][0], legend=True, edgecolor='black', cmap='Reds')
X_kde.plot(column='highway_bus_stop', ax=ax[2][1], legend=True, edgecolor='black', cmap='Greens')
X_kde.plot(column='highway_residential', ax=ax[2][2], legend=True, edgecolor='black', cmap='Blues')

ax[0][0].set_ylabel('PoI data', fontsize=18)
ax[1][0].set_ylabel('Quadrat', fontsize=18)
ax[2][0].set_ylabel('KDE', fontsize=18)

ax[2][0].set_xlabel('amenity_hospital', fontsize=18)
ax[2][1].set_xlabel('highway_bus_stop', fontsize=18)
ax[2][2].set_xlabel('highway_residential', fontsize=18)


sns.despine(trim=True)

fig.tight_layout()
fig.savefig('qualitative.pdf', dpi=200)<jupyter_output><empty_output><jupyter_text>Quantitive assessment Spatial heterogeneity (q) needs to use stratas on its calculation. We chose to get the city's neighborhoods as stratas.<jupyter_code>nhoods = api.get(city, admin_level=10)[['geometry', 'name']]
nhoods.plot();
s_grid = s_grid.join(gpd.sjoin(s_grid, nhoods)['name'].to_frame('strata'))
s_grid = gpd.GeoDataFrame(s_grid).drop_duplicates(subset=['geometry'])
s_grid.plot(column='strata');
scores_quadrat, scores_kde = {}, {}
for feature in X_kde.drop(columns=['geometry']).columns:
    scores_kde[feature] = {'q':gh.util.q_ongrid(X_kde[feature], s_grid, strata_col='strata'),
                           'I':gh.util.moran_i_ongrid(X_kde[feature], 
                                                      s_grid['geometry'].centroid.apply(lambda x:x.coords[0]), 
                                                      d_threshold=resolution/110.5)[0]}

    scores_quadrat[feature] = {'q':gh.util.q_ongrid(X_quadrat[feature], s_grid, strata_col='strata'),
                               'I':gh.util.moran_i_ongrid(X_quadrat[feature], 
                                                          s_grid['geometry'].centroid.apply(lambda x:x.coords[0]), 
                                                          d_threshold=resolution/110.5)[0]}
scores_kde = scores_kde = pd.DataFrame(scores_kde).T.reset_index()\
                    .melt(value_vars=['q','I'], id_vars=['index'])\
                    .assign(Method='KDE')
scores_quadrat = pd.DataFrame(scores_quadrat).T.reset_index()\
                    .melt(value_vars=['q','I'], id_vars=['index'])\
                    .assign(Method='Quadrat')

scores = pd.concat([scores_quadrat, scores_kde]).rename(columns={'index':'feature'})
fig, ax = plt.subplots(ncols=2, figsize=(12, 5))
sns.barplot(y='feature', x='value', hue='Method', data=scores.loc[scores['variable']=='q'], ax=ax[0], palette='viridis')
sns.barplot(y='feature', x='value', hue='Method', data=scores.loc[scores['variable']=='I'], ax=ax[1], palette='viridis')


ax[0].yaxis.set_tick_params(labelsize=14)
ax[0].xaxis.set_tick_params(labelsize=12)
ax[0].set_ylabel('')
ax[0].set_xlabel('Spatial heterogeneity\nq', fontsize=18)
ax[0].get_legend().remove()

ax[1].legend(fontsize=18)
ax[1].xaxis.set_tick_params(labelsize=12)
ax[1].set_yticks([]) 
ax[1].set_ylabel('', fontsize=0)
ax[1].set_xlabel('Spatial autocorrelation\nI', fontsize=18)
fig.tight_layout()

fig.savefig('quantitative.pdf', dpi=200)<jupyter_output><empty_output><jupyter_text>Demonstration and applications<jupyter_code>from sklearn.cluster import KMeans

features = ['amenity_hospital', 'highway_bus_stop', 'highway_residential']
X = X_kde[features]

kmeans = KMeans(3).fit(X)
X['Clusters'] = kmeans.labels_
X['Clusters'] = X['Clusters'].apply(lambda x: [f'C{i}' for i in range(5)][x])
fig = sns.pairplot(X, hue='Clusters', size=1.5, palette='Set1')
fig.savefig('cluster_histograms.png', dpi=500)
norm=plt.Normalize(-2,2)
cmap = matplotlib.colors.LinearSegmentedColormap.from_list("", ["#E53631","#377EB8","#50AF4A"])

fig, ax = plt.subplots()
X['geometry'] = X_kde['geometry']
gpd.GeoDataFrame(X).plot(column='Clusters', cmap=cmap, edgecolor='black', alpha=0.9, legend=True, ax=ax)
sns.despine(trim=True)
fig.savefig('clusters_map.png', dpi=200)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1943900848010437, 'fraction_numerical': 0.020743639921722113, 'mean_word_length': 4.924265842349304, 'pattern_counts': {'":': 0, '<': 36, '<?xml version=': 0, '>': 36, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '491741', 'license': 'MIT', 'n_tokens_mistral': 3031, 'n_tokens_neox': 2881, 'n_words': 352, 'path': 'examples/paper_experiment.ipynb', 'repo_name': 'adaj/geohunter', 'size': '7665'}