Une introduction immersive à l'analyse géospatiale

Jan 16 2023

Recherche d'appartements idéalement situés à Singapour En 2021, mon équipe et moi avons remporté un concours de visualisation de données avec plus de 400 équipes participantes. Notre soumission a exploré une nouvelle façon de rechercher des appartements, en se concentrant sur l'accessibilité des commodités comme critère de recherche principal.

Recherche d'appartements idéalement situés à Singapour

Immeubles d'appartements près du centre-ville de Singapour (Photo de Mike Enerio sur Unsplash)

En 2021, mon équipe et moi avons remporté un concours de visualisation de données avec plus de 400 équipes participantes. Notre soumission a exploré une nouvelle façon de rechercher des appartements, en se concentrant sur l'accessibilité des commodités comme critère de recherche principal.

Ce projet m'a initié à une gamme de techniques d'analyse géospatiale que je partagerai dans cet article. J'ai également partagé les ensembles de données et le code Python utilisés pour notre soumission sur ma page GitHub afin que vous puissiez suivre une introduction immersive à l'analyse géospatiale.

Table des matières : (cliquez pour passer à la section)
1. Jeux de données géospatiales
2. Géocodage
3. Systèmes de coordonnées de référence
4. Visualisation
5. Analyse
6. Rassembler le tout

Ensembles de données géospatiales

Les ensembles de données géospatiales intègrent des informations spatiales (où se trouve quelque chose) et des informations d'attribut (ce qu'est quelque chose) . Par exemple, la couche de trafic dans Google Maps repose sur un ensemble de données qui intègre l'emplacement et la vitesse du trafic le long de chaque route.

Instantané de la couche de trafic Google Maps pour Dallas, Texas

Les ensembles de données géospatiales se présentent sous divers formats. Dans cet article, je vais aborder la façon de travailler avec deux que j'ai fréquemment rencontrés : les fichiers de formes et les fichiers KML. Dans les deux cas, nous utiliserons la bibliothèque python Geopandas .

Fichiers de formes (.shp)

Le format Shapefile a été développé par Esri, un leader du marché des logiciels d'analyse géospatiale. Les fichiers de formes sont compatibles avec les logiciels Esri comme ArcGIS ainsi qu'avec d'autres plateformes comme Python et Tableau.

Lors de la recherche d'appartements à Singapour, la proximité des gares est une considération clé car c'est le principal mode de transport. Vous pouvez télécharger l'emplacement des gares sous forme de shapefile ici . Vous remarquerez que le jeu de données comprend un ensemble de fichiers (.shp, .shx, .dbf, etc.), comme c'est le cas pour les fichiers de formes. Assurez-vous que ces fichiers conservent le même nom et sont stockés dans le même dossier, sinon ils ne s'ouvriront pas correctement.

Les fichiers de formes peuvent être lus à l'aide de la fonction geopandas read_file .

import geopandas as gpd
filepath = "<insert path to folder with shapefile>/MRTLRTStnPtt.shp"
mrt = gpd.read_file(filepath)

      
                
                    
                    
       First 3 lines of mrt GeoDataframe

Fichiers Keyhole Markup Language (.kml)

Le format KML a été développé à l'origine pour être utilisé avec Google Earth, mais est depuis devenu une norme internationale. L'un de leurs principaux avantages est qu'ils peuvent stocker un mélange de différentes géométries (par exemple, des points et des polygones), alors que les fichiers de formes ne prennent en charge que le type de géométrie homogène.

KML utilise une structure basée sur des balises avec des éléments imbriqués. Par exemple, ce qui suit est un extrait d'un fichier KML de l'emplacement des supermarchés à Singapour ( téléchargez ici ). Chaque paire de balises <Placemark id=XXX>…</Placemark> relie les données d'un objet géographique, dans ce cas un supermarché. Les balises imbriquées sont utilisées pour étiqueter les différents champs de données pour cet objet.

Extrait du fichier KML des supermarchés à Singapour. Les surlignages indiquent les balises. Les couleurs différencient les niveaux d'imbrication : les balises bleues sont imbriquées dans les balises jaunes, elles-mêmes imbriquées dans les balises vertes.

Lors de la lecture de fichiers kml à l'aide de Geopandas, seules certaines de ces balises sont automatiquement analysées en colonnes.

import geopandas as gpd
gpd.io.file.fiona.drvsupport.supported_drivers['KML'] = 'rw' # activate KML driver
filepath = "<insert path to folder with kml file>/supermarkets-kml.kml"
sup_mkt = gpd.read_file(filepath,driver="KML")

      
                
                    
                    
       First 3 lines of sup_mkt GeoDataframe. ‘Description’ column requires further parsing

def parse_supmkt(string):
    tags = ['LIC_NAME', 'BLK_HOUSE', 'STR_NAME', 'UNIT_NO', 'POSTCODE']
    
    output = []
    
    for tag in tags:
        start = "<th>" + tag + "</th> <td>"
        end = "</td>"
        len_start = len(start)
        
        i_start = string.find(start)
        i_end = string.find(end,i_start+len_start)
        
        entry = string[(i_start+len_start):i_end]
        output.append(entry)
    
    return output

supmkt_1 = sup_mkt.copy()
supmkt_1['Lic_Name'] = supmkt_1['Description'].apply(lambda x: parse_supmkt(x)[0])
supmkt_1['Address'] = supmkt_1['Description'].apply(lambda x: parse_supmkt(x)[1] + " " + parse_supmkt(x)[2]
                                                  + ", #" + parse_supmkt(x)[3])
supmkt_1['Postcode'] = supmkt_1['Description'].apply(lambda x: parse_supmkt(x)[4])
supmkt_1.drop(columns=['Name','Description'],inplace=True)

      
                
                    
                    
       First 3 lines of supmkt_1 GeoDataframe. Licensee name, Address and Postcode have now been parsed.

Le géocodage est la conversion d'adresses en coordonnées géographiques (ou vice versa) . Les ensembles de données brutes décrivent parfois l'emplacement des points d'intérêt en fonction de leur adresse postale. Un exemple est l' ensemble de données HDB Property Information qui fournit des informations sur les immeubles d'appartements à Singapour.

Extrait de l'ensemble de données HDB Property Information. L'emplacement est décrit à l'aide de l'adresse municipale (boîte rouge).

Pour une analyse géospatiale précise telle que le calcul des distances entre les points, ces adresses doivent être converties en coordonnées géographiques. Cette conversion se fait généralement en transmettant les adresses postales à un service de géocodage en ligne à l'aide de leur API.

Géocodage avec l'API OneMap

Pour les adresses à Singapour, l' API de recherche OneMap fournie par la Singapore Land Authority est un service gratuit et précis.

Essentiellement, les coordonnées de chaque immeuble peuvent être obtenues en passant l'adresse de la rue (concaténation des colonnes 'blk_no' et 'street' de notre jeu de données) à la fonction ci-dessous :

import requests 

def getcoordinates(address):
    req = requests.get('https://developers.onemap.sg/commonapi/search?searchVal='+address+'&returnGeom=Y&getAddrDetails=Y&pageNum=1')
    resultsdict = eval(req.text) #evaluate as python expression
    #keep only those with results
    if len(resultsdict['results'])>0:
        return resultsdict['results'][0]['LATITUDE'], resultsdict['results'][0]['LONGITUDE'], resultsdict['results'][0]['POSTAL'] #i believe the '[0]' returns the top entry
    else:
        pass

from tqdm import tqdm # to view progress bar 

coordinateslist= []
failed_count = 0
addresslist = resi_addresses['addr'] # 'addr' is a concatenation of 'blk_no' and 'street'
print(f"Total no. of addresses: {len(addresslist)}")
for address in tqdm(addresslist, desc="Geocoding..."): # update progress bar at each iteration
    try: #try executing the following
        coordinates = getcoordinates(address)
        if len(coordinates)>0:
            coordinateslist.append(coordinates)
            
    except: #if you receive an error, do this          
        failed_count += 1
        print(f"No. of addresses with no coordinates (cummulative): {failed_count}")
        coordinateslist.append(None)

df_coordinates = pd.DataFrame(coordinateslist)
df_combined = resi_addresses.join(df_coordinates)
df_combined  = df_combined .rename(columns={0:'Latitude', 1:'Longitude', 2:'Postcode'}).drop(columns = ['addr'])

      
                
                    
                    
       3 lines from df_combined. ‘blk_no’ and ‘street’ are from the original HDB Property Info dataset. ‘Latitude’, ‘Longitude’ and ‘Postcode’ were retrieved through geocoding.

De nombreuses autres API de géocodage sont disponibles. La bibliothèque Python Geopy vous permet d'implémenter plusieurs de ces API en utilisant une syntaxe standardisée, au lieu d'apprendre la syntaxe unique de chaque API.

Nominatim est un service de géocodage gratuit disponible via Geopy. Il est plus polyvalent que OneMap en raison de sa couverture mondiale, mais j'ai choisi d'utiliser OneMap pour ce projet car il était plus précis pour les adresses à Singapour.

Systèmes de référence de coordonnées

Jusqu'à présent, dans cet article, j'ai fait référence aux coordonnées géographiques de manière assez vague. Il est important de noter qu'il existe de nombreux systèmes de référence de coordonnées (CRS) différents et que les valeurs de coordonnées doivent être interprétées en fonction du CRS sur lequel elles sont basées . Différents CRS ont différentes dispositions de grille xy, unités, données/origine, etc.

Une illustration simple de la façon dont les mêmes valeurs de coordonnées (1,1) définies sur la base de différents systèmes de coordonnées correspondent en réalité à différents emplacements.

Il existe deux grandes classes de CRS : CRS géographique et projeté . Le tableau ci-dessous les compare.

Identification des CRS existants

Les ensembles de données géospatiales comme les fichiers de formes et les fichiers kml sont généralement intégrés aux informations CRS. Une fois importé en tant que GeoDataframe, leur SCR peut être lu facilement :

Lire le CRS de l'ensemble de données des supermarchés avec lequel nous avons travaillé précédemment

Affectation de CRS

Parfois, les ensembles de données géospatiales ne sont pas fournis avec un CRS intégré. Par exemple, lorsque nous avons géocodé le jeu de données HDB Property Info, nous n'avons pas intégré les informations CRS dans le dataframe. Pour ce faire, nous devons d'abord créer une colonne 'géométrie', pour stocker les informations spatiales de chaque objet en tant que classe de géométrie Shapely.

import pandas as pd
from shapely.geometry import Point

filepath = "<insert filepath for geocoded HDB Property Info Dataset>" 
resi_blks = pd.read_csv(filepath, usecols=['blk_no','street','Latitude','Longitude','Postcode'])

# Represent each apartment block as a Shapely geometry (point)
resi_blks['geometry'] = resi_blks.apply(lambda row: Point(row.Longitude, row.Latitude), axis=1)

      
                
                    
                    
       First 3 lines of resi_blks dataframe, including the newly created ‘geometry’ column

import geopandas as gpd

resi_blks2 = gpd.GeoDataFrame(resi_blks,crs="EPSG:4326") #convert to GeoDataframe
resi_blks2.to_file("resi_blks2.shp") # save as shapefile

La reprojection fait référence à la conversion des coordonnées géographiques d'un CRS à un autre. Dans l'exemple ci-dessous, nous reprojetons nos informations de propriété HDB GeoDataframe d' EPSG 4326 (un CRS géographique) à EPSG 32648 (un CRS projeté adapté à Singapour) pour faciliter une mesure de distance plus précise.

resi_blks2 = resi_blks2.to_crs("EPSG:32648")

Pour une visualisation étendue, la création de tableaux de bord, etc., vous souhaiterez peut-être utiliser un logiciel plus approprié comme Tableau ou ArcGIS. Cependant, ceux-ci ne sont pas gratuits. En outre, il est souvent utile de visualiser vos données dans l'environnement python afin de pouvoir vérifier que les étapes de géotraitement effectuées renvoient des résultats raisonnables . Nous en verrons des exemples dans la section suivante.

Geopandas a une méthode de tracé qui peut être directement appliquée à n'importe quel GeoDataframe. Cependant, cela n'est généralement pas significatif car il n'y a pas de carte de base pour donner un contexte spatial.

Emplacement des gares tracé à l'aide de geopandas.plot().

Cartographie interactive avec Folium

La bibliothèque Folium vous permet de tracer vos données sur une carte de base, et permet une cartographie interactive (par exemple, vous pouvez zoomer, déplacer) avec une carte de base . Il existe de nombreuses fonctions pour prendre en charge différentes personnalisations. Je ne couvrirai que l'essentiel dans cet article.

Commencez par créer une carte de base à l'aide de la fonction Map de Folium . La sélection d'un emplacement approprié pour centrer votre carte et le niveau de zoom de départ nécessitera quelques essais et erreurs. Dans cet exemple, ils ont été ajustés pour montrer Singapour.

import folium

map1 = folium.Map(location=([1.3521,103.8198]), # [latitude, longitude]. CRS: EPSG 4326 CRS
                  zoom_start = 11.5) #initial zoom level
map1 # display

      
                
                    
                    
       map1 (base map only). By default, OpenStreetMap map tiles are used. Other options are available.

for _,row in mrt.iterrows(): 
    row['location'] = list(reversed((row['geometry'].x,row['geometry'].y))) # folium requires location as [latitude, longitude]
    folium.map.Marker(location=row['location'],            # marker location
                     icon=folium.Icon(color='red',         # custom marker style
                                     icon_color='black',
                                     icon='train',
                                     prefix='fa'),         # use this prefix for fontawesome icons
                     popup=row['STN_NAME']).add_to(map1)   #info to display in pop-up
map1 # display

      
                
                    
                    
       map1 with markers at train stations. Top left: full map; bottom right: zoomed in view showing a pop-up.

Analyse

La bibliothèque Geopandas dispose d'un large éventail de fonctions que vous pouvez utiliser pour différents types d'analyse. Pour vous donner une idée de ce qui est possible, voici quelques exemples :

Génération de tampons

La méthode du tampon Geopandas génère une nouvelle géométrie qui englobe tous les points situés à une distance spécifiée de la géométrie d'entrée. Pour ce projet, des zones tampons ont été utilisées pour définir la zone dans un certain rayon autour de chaque immeuble en tant que précurseur d'une analyse plus approfondie.

Rappelez-vous l'ensemble de données géocodé HDB Property Info que nous avons vu précédemment :

3 premières lignes de resi_blks2 GeoDataframe.

Utilisez le code ci-dessous pour générer des tampons de rayons 400 m — 2 km autour de chaque immeuble .

# use projected CRS as buffer generation involves distance measurement
resi_blks2 = resi_blks2.to_crs("EPSG:32648")

buffer_distances = [400,800,1200,1600,2000]

resi_blks3 = resi_blks2.copy()

for i in range(len(buffer_distances)):
    resi_blks3['buffer_' + str(buffer_distances[i])] = resi_blks3['geometry'].buffer(buffer_distances[i])

      
                
                    
                    
       First 3 lines of resi_blks3

test_buffers = resi_blks3.loc[[0,1],:] #extract first 2 blocks

# Convert crs to epsg 4326 to match what Folium supports
test_buffers = test_buffers.to_crs("EPSG:4326")
for i in range(len(buffer_distances)):
    test_buffers = test_buffers.set_geometry('buffer_' + str(buffer_distances[i])).to_crs("EPSG:4326")

map1 = folium.Map(location=([1.3521,103.8198]),zoom_start=12) #folium basemap

for _,row in test_buffers.iterrows():
    #add markers for blk locations
    row['location'] = (row['geometry'].x,row['geometry'].y)
    folium.map.Marker(location=list(reversed(row['location'])),
                     icon=folium.Icon(color='red',
                                     icon_color='black',
                                     icon='home'),
                     popup=row['blk_no']+" "+row['street']).add_to(map1)
    
    #add buffers
    for i in range(len(buffer_distances)):
        buffer_gs = gpd.GeoSeries(row['buffer_' + str(buffer_distances[i])])
        buffer_json = buffer_gs.to_json()
        folium.GeoJson(data=buffer_json).add_to(map1)

map1 # display map

      
                
                    
                    
       map1 showing 400m — 2km buffers generated around two apartment blocks

La méthode Geopandas within détermine si une géométrie est contenue dans une autre. Nous pouvons l'utiliser pour compter les commodités contenues dans les zones tampons de 400 m (5 à 10 minutes de marche) autour de chaque immeuble.

À titre d'exemple, nous utiliserons l'ensemble de données des supermarchés que nous avons vu précédemment dans cet article :

3 premières lignes de supmkt GeoDataframe

Utilisez le code ci-dessous pour compter le nombre de supermarchés dans un rayon de 400 m autour de chaque appartement.

def count_amenity(amenity_gdf, poly_buffer):
    """
    Compute number of amenities within the specified buffer polygon
    amenity_gdf = geodataframe for amenity type of interest
    polygon_buffer = buffer polygon to search within
    """
    
    return (len(amenity_gdf.loc[amenity_gdf.within(poly_buffer)]))

amenities_1 = resi_blks3.copy()
supmkt = supmkt.to_crs(amenities_1.crs) # ensure same crs for meaningful overlay
amenities_1['n_supmkt_400m'] = amenities_1['buffer_400'].apply(lambda x: count_amenity(supmkt, x))

      
                
                    
                    
       First 3 lines of amenities_1

La carte ci-dessous montre l'emplacement des espaces verts à travers Singapour ( téléchargez ici ).

Localisation des espaces verts à Singapour

La méthode de superposition Geopandas génère de nouvelles géométries en fonction de la façon dont deux ensembles de géométries d'entrée se chevauchent (ou ne se chevauchent pas).

Illustration de différentes options de superposition à l'aide de Geopandas. (Source : Documentation Géopandas)

Pour ce projet, nous pouvons utiliser l' option d'intersection pour identifier les zones vertes situées dans les zones tampons de 400 m (5 à 10 minutes de marche) autour de chaque immeuble. On peut alors utiliser la méthode des surfaces de Geopandas pour calculer la surface de ces espaces verts identifiés. L'image ci-dessous illustre cela pour un immeuble.

Espaces verts (polygones verts) à moins de 400 m (cercle bleu) d'un immeuble (marqueur rouge)

Utilisez le code suivant pour exécuter le calcul pour tous les immeubles :

def total_green_area(poly_buffer): 
    """Sum up total green area (in ha) within specified polygon buffer (poly_buffer)"""
    
    #convert poly_buffer to geodataframe
    polybuffer_gdf = gpd.GeoDataFrame({'geometry':[poly_buffer]},crs=greenery.crs)
    
    #find the intersection polygons between poly_buffer and greenery geodataframes
    green_spaces = gpd.overlay(polybuffer_gdf,greenery,how='intersection')
    
    #calculate area of each green space
    green_spaces['area'] = green_spaces['geometry'].area
    
    total_area = green_spaces['area'].sum()
    total_area = total_area/10000 #convert m2 to ha
    
    return (total_area)

amenities_2 = amenities_1.copy()
amenities_2['grnarea_400'] = amenities_2['buffer_400'].apply(lambda x: total_green_area(x))

      
                
                    
                    
       First 3 lines of amenities_2

La méthode de distance Geopandas calcule la distance entre deux géométries d'entrée. Nous pouvons l'utiliser pour calculer la distance entre les immeubles et chaque gare à partir de l'ensemble de données que nous avons vu précédemment :

3 premières lignes de mrt GeoDataframe

Nous pourrions calculer la distance entre chaque appartement et chaque gare et identifier la plus proche, mais ce serait très lent. Au lieu de cela, nous limitons d'abord notre recherche à la zone tampon de 400 m autour de chaque appartement. S'il n'y a pas de stations dans cette zone, nous élargissons notre recherche à la zone tampon de 800 m, et ainsi de suite progressivement jusqu'à un maximum de 2 km.

def dist_to_amenity(amenity_gdf, row):
    """
    Compute distance to nearest amenity
    amenity_gdf = geodataframe for amenity type of interest
    row = row of data with various buffer polygons to search within
    """
    #list amenities within 400m
    mybuffer = row.buffer_400
    list_amen = amenity_gdf.loc[amenity_gdf.within(mybuffer)]
    
    #if none within 400m, widen buffer to 800m and so on
    if len(list_amen)==0:
        mybuffer = row.buffer_800
        list_amen = amenity_gdf.loc[amenity_gdf.within(mybuffer)]
        
        if len(list_amen)==0:
            mybuffer = row.buffer_1200
            list_amen = amenity_gdf.loc[amenity_gdf.within(mybuffer)]
        
            if len(list_amen)==0:
                mybuffer = row.buffer_1600
                list_amen = amenity_gdf.loc[amenity_gdf.within(mybuffer)]
                
                if len(list_amen)==0:
                    mybuffer = row.buffer_2000
                    list_amen = amenity_gdf.loc[amenity_gdf.within(mybuffer)]
    
    #if no amenity within 2km, return null for distance
    if len(list_amen)==0:
        return(np.nan)
    else:
        # compute distance from each amenity to apartment block
        list_amen['distance'] = list_amen['geometry'].distance(row.geometry)
        return(min(list_amen['distance'])) # return minimum distance

amenities_3 = amenities_2.copy()
mrt = mrt.to_crs(amenities_3.crs) # ensure same CRS for meaningful overlay
amenities_3['dist_mrt']=amenities_3.apply(lambda row: dist_to_amenity(mrt, row), axis=1)

      
                
                    
                    
       First 3 lines of amenities_3

L'objectif de ce projet était d'explorer une nouvelle façon de rechercher des appartements à Singapour, en se concentrant sur l'accessibilité des commodités comme critère de recherche principal. Dans cet article, j'ai montré comment cartographier l'emplacement des appartements à Singapour par géocodage et comment extraire l'emplacement des commodités à partir d'ensembles de données géospatiales. En utilisant les techniques d'analyse discutées dans la dernière section, nous pouvons ensuite quantifier et comparer l'accessibilité des différentes commodités de chaque appartement.

Pour afficher le code Python complet et le classeur Tableau utilisé, visitez ma page GitHub .

Les références

Introduction aux données vectorielles spatiales en Python (par Earth Data Science)
Géocodage des coordonnées de Singapour : API OneMap (par Shawn Tham, Towards Data Science)