Zones isochrones
Rendez-vous sur https://openrouteservice.org/ et créez-vous un compte afin de récupérer une clé d'API (le token).
Il s'agit du service associé au plugin ORS Tools, que vous pouvez aussi installer afin de faire quelques tests.
Limitations
Mais le plugin est bridé à 500 zones isochrones par jour et surtout, à des batchs de 5 zones par traitement et 60 zones par minute. Or nous avons 300 et quelques zones isochrones à créer...
Nous allons essayer d'outrepasser ces limites en créant nos zones isochrones dans une boucle Python.
Construction d'un code de création d'une seule zone isochrone
Rendez-vous ici (une fois connecté au service) pour paramétrer la création d'une zone isochrone de test, la visualiser et en extraire son code Python. Utilisez les coordonnées du magasin de Cergy (2.0495375241349554, 49.05473127094994) afin de vérifier que vos paramètres sont bons.
import requests body = {"locations":[[2.0495375241349554,49.05473127094994]],"range":[900],"attributes":["area","reachfactor","total_pop"],"range_type":"time"} headers = { 'Accept': 'application/json, application/geo+json, application/gpx+xml, img/png; charset=utf-8', 'Authorization': 'YOUR-API-KEY', 'Content-Type': 'application/json; charset=utf-8' } call = requests.post('https://api.openrouteservice.org/v2/isochrones/driving-car', json=body, headers=headers) print(call.status_code, call.reason) print(call.text)
Exécutez ce code dans la console Python de QGIS. Le dernier print vous renvoie un GEOJSON (sous forme textuelle pour l'instant).
Vous avez déjà visualisé votre zone isochrone dans la cartographie d'ORS, mais vérifiez que le GEOJSON généré par le code dans QGIS est lui-même valide (en le copiant-collant dans l'outil geojson.io, ou en l'enregistrant puis en l'ouvrant dans QGIS).
Observez les coordonnées que nous utilisons dans le code : elles sont en 4326 (WGS84).
Affichage sous QGIS
Si tout va bien, commentez les print et ajoutez à la suite de votre code les lignes suivantes pour afficher le GEOJSON dans QGIS (en modifiant les chemins) :
with open('C:/Users/Georges/Downloads/myfile.geoson', 'w', encoding='utf-8') as outfile: outfile.write(call.text) iface.addVectorLayer('C:/Users/Georges/Downloads/myfile.geoson', 'Isochrone', 'ogr')
Exécutez l'ensemble du code dans la console Python de QGIS. La zone isochrone créé avec Python apparaît.
Itération
Bien, nous savons créer une zone isochrone dans QGIS avec Python, mais il nous en faut 300 et quelques !
Nous allons donc devoir itérer à travers la couche de points pour créer les zones isochrones correspondantes à partir de leur coordonnées en 4326. Nous ajoutons donc une alerte sur la projection de la couche (if...).
Pour nos tests, nous ne travaillerons que sur 3 points (limit = 3).
import requests mylayer = QgsProject.instance().mapLayersByName("decathlon_france")[0] layerProj = mylayer.crs() limit = 3 i = 0 for feat in mylayer.getFeatures(): id_mag = feat['id'] geom = feat.geometry() i += 1 if i > limit: break print('\nMagasin ' + str(id_mag) + ' (itération ' + str(i) + ')') print('X : ' + str(geom.asPoint().x())) print('Y : ' + str(geom.asPoint().y())) if '4326' not in str(layerProj): print('\nATTENTION : projection différente de 4326 :') else: print('\nOK : la couche est en 4326 :') print(str(layerProj) + '\n')
Reprojection
Si besoin est, nous allons utilisé ce code ci-dessous pour :
- reprojeter notre couche
- l'enregistrer dans un autre shape suffixé "4326"
- fermer la couche d'origine
- ouvrir la nouvelle couche en la nommant de l'exacte même façon que la couche originale. Ainsi nos autres codes fonctionneront de l'exacte même façon.
project = QgsProject.instance() mylayer = QgsProject.instance().mapLayersByName("decathlon_france")[0] mylayer4326 = 'C:/Users/georg/Downloads/decathlon_france4326.shp' processing.run("native:reprojectlayer", { 'INPUT':mylayer, 'TARGET_CRS':QgsCoordinateReferenceSystem('EPSG:4326'), 'CONVERT_CURVED_GEOMETRIES':False, 'OPERATION':'+proj=pipeline +step +inv +proj=lcc +lat_0=46.5 +lon_0=3 +lat_1=49 +lat_2=44 +x_0=700000 +y_0=6600000 +ellps=GRS80 +step +proj=unitconvert +xy_in=rad +xy_out=deg', 'OUTPUT':mylayer4326 }) project.addMapLayer(mylayer) project.removeAllMapLayers() mylayerOpen = QgsVectorLayer(mylayer4326, 'decathlon_france', 'ogr') project.addMapLayer(mylayerOpen)
Marquer les magasins traités
Bien, maintenant il s'agit de marquer les magasins dont la zone a déjà été créée. Car pour contourner les limitations de l'API, et pouvoir relancer notre code à souhait, il nous faudra limiter la boucle à quelques dizaines de points, et marquer les magasins dont la zone a déjà été créée. Ainsi nous n'aurons qu'à ré-exécuter le code plusieurs fois d'affilé pour générer nos 300 et quelques zones isochrones. Attention à la limite de 500 zones par jour, qui elle ne pourra pas être dépassée.
Ajoutez un champ iso dans votre couche de points de test, en integer, à 0. Nous le faisons passer à 1 à chaque zone isochrone créée, et ajoutons une condition pour ne passer que sur les magasins dont la zone n'a pas encore été créée :
with edit(mylayer): for feat in mylayer.getFeatures(): if feat['iso'] == 0: id_mag = feat['id'] geom = feat.geometry() i += 1 if i > limit: break feat['iso'] = 1 mylayer.updateFeature(feat) print(id_mag, geom.asPoint().x(), geom.asPoint().y(), i)
Maintenant, à la suite de la boucle, la zone isochrone, où nous remplaçons les coordonnées en dur par nos géométries, avec l'id des magasins en nom de fichier :
body = {"locations":[[geom.asPoint().x(), geom.asPoint().y()]],"range":[900],"attributes":["area","reachfactor","total_pop"],"range_type":"time"} headers = { 'Accept': 'application/json, application/geo+json, application/gpx+xml, img/png; charset=utf-8', 'Authorization': 'YOUR-API-KEY', 'Content-Type': 'application/json; charset=utf-8' } call = requests.post('https://api.openrouteservice.org/v2/isochrones/driving-car', json=body, headers=headers) #print(call.status_code, call.reason) #print(call.text) with open('C:/Users/Georges/Downloads/CoursQGIS/data/isochrones/'+id_mag+'.geoson', 'w', encoding='utf-8') as outfile: outfile.write(call.text) iface.addVectorLayer('C:/Users/Georges/Downloads/CoursQGIS/data/isochrones/'+id_mag+'.geoson', id_mag, 'ogr')
Mais nous n'avons aucun attribut pour identifier les zones isochrones par magasin. Et nous en aurons besoin par la suite.
Nous ajoutons donc un champ dans le fichier GEOJSON, que nous remplissons avec l'id des magasins.
my_geojson = QgsProject.instance().mapLayersByName(id_mag)[0] pr = my_geojson.dataProvider() pr.addAttributes([QgsField("id_mag", QVariant.String)]) my_geojson.updateFields() with edit(my_geojson): for feature in my_geojson.getFeatures(): feature['id_mag'] = id_mag my_geojson.updateFeature(feature)
OK. Testez le code sur votre extraction de 10 magasins, si cela fonctionne, vous n'avez plus qu'à remplacer la variable mylayer et à augmenter la limit (n'oubliez pas de supprimer les zones précédemment apparues, à la fois dans QGIS mais aussi les fichiers GEOJSON dans votre machine, et de créer le champ iso dans la couche de vos magasins).
Code complet
import requests mylayer = QgsProject.instance().mapLayersByName("decathlon_france")[0] limit = 20 i = 0 with edit(mylayer): for feat in mylayer.getFeatures(): if feat['iso'] == 0: id_mag = feat['id'] geom = feat.geometry() i += 1 if i > limit: break feat['iso'] = 1 mylayer.updateFeature(feat) # print(id_mag, geom.asPoint().x(), geom.asPoint().y(), i) # ZONE ISOCHRONE body = {"locations":[[geom.asPoint().x(), geom.asPoint().y()]],"range":[900],"attributes":["area","reachfactor","total_pop"],"range_type":"time"} headers = { 'Accept': 'application/json, application/geo+json, application/gpx+xml, img/png; charset=utf-8', 'Authorization': 'YOUR-API-KEY', 'Content-Type': 'application/json; charset=utf-8' } call = requests.post('https://api.openrouteservice.org/v2/isochrones/driving-car', json=body, headers=headers) #print(call.status_code, call.reason) #print(call.text) with open('C:/Users/Georges/Downloads/CoursQGIS/data/isochrones/'+id_mag+'.geoson', 'w', encoding='utf-8') as outfile: outfile.write(call.text) iface.addVectorLayer('C:/Users/Georges/Downloads/CoursQGIS/data/isochrones/'+id_mag+'.geoson', id_mag, 'ogr') # AJOUTER ID MAGASIN my_geojson = QgsProject.instance().mapLayersByName(id_mag)[0] pr = my_geojson.dataProvider() pr.addAttributes([QgsField("id_mag", QVariant.String)]) my_geojson.updateFields() with edit(my_geojson): for feature in my_geojson.getFeatures(): feature['id_mag'] = id_mag my_geojson.updateFeature(feature)
Hop ! Vous obtenez 300 et quelques zones isochrones à 15 minutes en voiture pour atteindre un magasin Décathlon ! Merci aux contributeurs OpenStreetMap et à l'université d'Heidelberg !
 isochrones15m.zip | [Les zones isochrones fusionnées] | 979 Ko |
| isochrones_ORS.zip | [Les fichiers GEOJSON des zones isochrones] | 786 Ko |
  Export OSM des magasins de sports français le 27 décembre 2020 | [key=shop ; value=sport ; in=France] | 0 Ko |