Principes de base de Python x SIG (partie 2)

Le premier est ici La partie 3 est ici Matériel pédagogique de l'Université d'Helsinki Nous résumerons les réponses et les suppléments pour la Semaine3-Semaine4.

Week3

3-1 Géocodage du centre commercial

L'objectif est de trouver le nombre d'habitants vivant à 1,5 km autour d'un grand centre commercial à Helsinki, c'est-à-dire la population de la zone commerciale. Vous devrez chercher vous-même sur le net pour connaître l'adresse du centre commercial. Présentation du géocodage pour convertir les adresses en coordonnées. Au Japon, Service de mise en correspondance d'adresses CSV fourni par l'Université de Tokyo > Est célèbre. De plus, des opérations familières à QGIS telles que la mise en mémoire tampon et le couplage spatial apparaîtront également.

import geopandas as gpd
import pandas as pd
import matplotlib.pyplot as plt
import requests
import geojson
from shapely.geometry import Polygon, LineString, Point
from pyproj import CRS
import os

#Lecture des données
data = pd.read_table('shopping_centers.txt', sep=';', header=None)
data.index.name = 'id'
data.columns=['name', 'addr']
#Géocodage
geo = geocode(data['addr'], provider='nominatim', user_agent='autogis_xx', timeout=4)
#Combiner des données
geo = geo.to_crs(CRS.from_epsg(3879))
geodata = geo.join(data)

#Tampon
geodata['buffer']=None
geodata['buffer'] = geodata['geometry'].buffer(distance=1500)
geodata['geometry'] = geodata['buffer']

#Acquisition de données de grille de population
url = 'https://kartta.hsy.fi/geoserver/wfs'
params = dict(service='WFS',version='2.0.0',request='GetFeature',
              typeName='asuminen_ja_maankaytto:Vaestotietoruudukko_2018',outputFormat='json')
r = requests.get(url, params=params)
pop = gpd.GeoDataFrame.from_features(geojson.loads(r.content))

#Conversion du système de coordonnées
pop = pop[['geometry', 'asukkaita']]
pop.crs = CRS.from_epsg(3879).to_wkt()
geodata = geodata.to_crs(pop.crs)

#Couplage spatial
join = gpd.sjoin(geodata, pop, how="inner", op="intersects")
#Calcul de la population de la zone commerciale
grouped = join.groupby('name')
for key, group in grouped:
    print('store: ', key,"\n", 'population:', sum(group['asukkaita']))

3-2 Centre commercial le plus proche

Trouvez le centre commercial le plus proche de votre domicile et de votre lieu de travail. Sauf si vous habitez en Finlande, indiquez un endroit approprié à Helsinki comme adresse de base. (L'importation de la bibliothèque est omise ci-dessous.)

#Lire les données
home = pd.read_table('activity_locations.txt', sep=';', header=None)
home.index.name='id'
home.columns = ['name', 'addr']
shop = pd.read_table('shopping_centers.txt', sep=';', header=None)
shop.index.name = 'id'
shop.columns=['name', 'addr']
#Géocodage
geo_home = geocode(home['addr'], provider='nominatim', user_agent='autogis_xx', timeout=4)
geo_shop = geocode(shop['addr'], provider='nominatim', user_agent='autogis_xx', timeout=4)
#Recherche du magasin le plus proche
destinations = MultiPoint(list(geo_shop['geometry']))
for home in geo_home['geometry']:
    nearest_geoms = nearest_points(home, destinations)
    print(nearest_geoms[1])

Week4

4-1 Visualisation des données d'accessibilité

Visualisez l'accessibilité en combinant les données de temps de trajet et les données du réseau de métro.

#Lire les données
grid = gpd.read_file('data/MetropAccess_YKR_grid_EurefFIN.shp')
data = pd.read_table('data/TravelTimes_to_5944003_Itis.txt', sep=';')
data = data[["pt_r_t", "car_r_t", "from_id", "to_id"]]
#Jointure de données
data_geo = grid.merge(data, left_on='YKR_ID', right_on='from_id')
#Données invalides(-1)Exclusion
import numpy as np
data_geo = data_geo.replace(-1, np.nan)
data_geo = data_geo.dropna()
#Hiérarchisation des données
import mapclassify
bins = [15, 30, 45, 60, 75, 90, 105, 120, 135, 150, 165, 180]
classifier = mapclassify.UserDefined.make(bins = bins)
data_geo['pt_r_t_cl'] = data_geo[['pt_r_t']].apply(classifier)
data_geo['car_r_t_cl'] = data_geo[['car_r_t']].apply(classifier)

#Visualisation
fig = plt.figure(figsize=(10,10))
ax1 = fig.add_subplot(1, 2, 1)  #Transport public
data_geo.plot(ax=ax1, column='pt_r_t_cl')
ax1.set_title("Itis-Travel times by PT")
ax2 = fig.add_subplot(1, 2, 2)  #Voiture privée
data_geo.plot(ax=ax2, column='car_r_t_cl')
ax2.set_title("Itis-Travel times by Car")
plt.tight_layout()
plt.show()
fig.savefig('itis_accessibility.png')

Il sera affiché comme ceci.

4-2 Zone d'influence du centre commercial

Sur la base des données d'accessibilité obtenues en 4-1, nous visons à visualiser la zone de puissance en trouvant le centre commercial le plus proche sur chaque réseau.

#Lire les données
filepaths = glob.glob('data/TravelTimes*.txt')
for path in filepaths:
    data = pd.read_table(path, sep=';')
    data = data[['from_id', 'pt_r_t']]
    data = data.rename(columns={'from_id':'YKR_ID'})
    #Nom de colonne'pt_r_t_{store}'changer en
    newname = path.replace('data/TravelTimes_to_', '')
    newname = newname.replace('.txt', '')
    newname = re.sub('\d{7}_', '', newname)
    data = data.rename(columns={'pt_r_t':'pt_r_t_'+newname})
    grid = grid.merge(data) #Combiner des données
grid = gpd.read_file('data/MetropAccess_YKR_grid_EurefFIN.shp')

#Données invalides(-1)Exclusion
import numpy as np
grid = grid.replace(-1, np.nan)
grid = grid.dropna()
#Distance la plus courte au centre commercial de chaque grille ・ Nom du centre commercial
grid['min_t'] = None
grid['dominant_service'] = None
columns = ['pt_r_t_Ruoholahti', 'pt_r_t_Myyrmanni','pt_r_t_Itis', 'pt_r_t_Jumbo', 'pt_r_t_IsoOmena', 'pt_r_t_Dixi','pt_r_t_Forum']
mini = lambda row:row[columns].min()
idx = lambda row:row[columns].astype(float).idxmin()
grid['min_t'] = grid.apply(mini, axis=1)
grid['dominant_service'] = grid.apply(idx, axis=1)

4-3 Population des centres commerciaux

C'est une chose solide d'agréger les données de la grille en utilisant dissoudre et de trouver la population dans la sphère à l'intersection. La partie qui chevauche 4-2 est omise.

#4-Passez à l'étape 2,Créer des données pop, y compris des sphères

#Dissoudre et intersecter
dissolved = grid.dissolve(by = 'dominant_service')
pop = pop[['geometry', 'asukkaita']] #Limité aux données nécessaires
intersection = gpd.overlay(grid, pop, how='intersection')

#Regrouper par sphère pour trouver la population de la sphère
grouped = intersection.groupby('dominant_service')
for key, group in grouped:
    print(key, ':', sum(group['asukkaita']))