Pour un projet récent, nous avons eu besoin de faire des calculs géographiques. Petit problème, ce projet est en Plone, et n'a donc pas de base SQL. Il est toujours possible d'en ajouter une mais les calculs étant simples, installer et maintenir un PostGIS nous paraissait trop lourd. Nous avons donc cherché d'autres solutions.

Nous avons donc trouvé deux librairies. Une en javascript (Turf), et une en python (Shapely). Pour des raisons de performance et de précision (un des calculs étant plus précis avec un fichier plus gros, nous y reviendrons plus tard), nous avons choisi d'utiliser Shapely.

Nous avions trois besoins différents :

• Calculer la taille d'une zone
• Savoir si cette zone était dans un cercle prédéfini
• Savoir dans quel pays se trouve cette zone

La bibliothèque Shapely utilise le même moteur que PostGIS, elle est donc très performante et contient toutes les fonctionnalités dont nous aurons besoin. Elle ne permet pas par contre de gérer nativement les systèmes de projection. Nous aurons donc besoin d'adapter les zones à nos besoins.

Calculer la taille d'une zone

Nous devons tout d'abord charger nos données. Nous stockons les zones en GeoJSON, en utilisant le système de projection WGS 84 (le classique latitude/longitude).

Il existe plusieurs manières de charger des données dans Shapely. Shapely supportant la Python Geo Interface, nous pouvons facilement charger notre zone en utilisant directement le GeoJSON :

from shapely.geometry import shape

# geojson doit être un dictionnaire python
# Nous pouvons transformer une chaine en utilisant json.loads()
zone = shape(geojson)

La première manière de récupèrer la taille de la zone serait d'utiliser l'attribut area. Mais cela ne nous donnerait pas ce que nous voulons. En effet, rappelez-vous, Shapely ne gère pas les projections. Nous aurions donc l'aire calculée bêtement à partir des latitudes/longitudes. Nous devons donc transformer nos coordonnées en mètres avant.

Nous allons pour cela utiliser la bibliothèque pyproj ainsi que la méthode transform de Shapely.

import pyproj    
import shapely
from functools import partial

zone_metres = shapely.ops.transform(
    partial(
        pyproj.transform,
        pyproj.Proj(init='EPSG:4326'),  # EPSG:4326 est WGS 84
        pyproj.Proj(
            proj='aea',
            lat1=zone.bounds[1],
            lat2=zone.bounds[3]
        )
    ),
    geom
)

Nous pouvons maintenant calculer la taille de notre zone :

taille_zone = zone_metres.area  # en m²

Calculer si notre zone est dans un cercle ou non

Il faut tout d'abord charger notre cercle. Comme c'est une forme simple (un polygone en fait), nous pouvons la stocker directement en python.

# Remplacer ... par les coordonnées de votre polygon
cercle = shape({"type": "Polygon", "coordinates": [...]})

Il ne nous reste plus qu'à faire le calcul. Nous avons plusieurs possibilités en fonction du résultat recherché.

# Vrai si la zone est entièrement dans le cercle
zone.within(cercle)

# Vrai si la zone touche le cercle
zone.interesects(cercle)

# Vrai si le centroïde de la zone est dans le cercle
zone.centroid.within(cercle)

Calculer dans quel pays se trouve notre zone

Le principe est exactement le même que pour le cercle, sauf que nous allons le faire avec tous les pays pour trouver le bon. Commençons par charger les pays. Nous pouvons par exemple récupérer les limites des pays en requêtant Open Street Map. Heureusement, quelqu'un l'a déjà fait pour nous !

Nous allons prendre les pays avec les frontières maritimes afin de trouver les zones un peu limite que l'on pourrait avoir. Au niveau de la taille, plus le GeoJSON est précis, plus il est lourd évidemment. Mais il prendra également plus de temps à charger. Trouver la bonne résolution dépendra donc de votre besoin.

Nous pouvons donc charger ce GeoJSON. Il est préférable de ne le faire qu'une seule fois, les données ne changeant pas d'une fois à l'autre. Deux petites choses à noter :

• Nous stockerons les géométries dans un dictionnaire avec le code à trois chiffre en clé (`A3` dans les propriétés de notre GeoJSON). Nous ne pourrons pas sinon associer la géométrie à un pays.
• Nous vérifions que la géométrie soit valide (pas de lignes qui se croisent, etc). Sinon, nous essayons de la corriger avec `.buffer(0)`, ce qui permet de redéfinir les contours de notre géométrie en éliminant certains petits défauts au passage.

with open(geojson_path) as geojson_file:
    geojson = json.loads(geojson_file.read())
liste_pays = {}
for pays in geojson.get('features', []):
    if 'geometry' in pays and 'A3' in pays.get('properties', {}):
        geom = shape(pays['geometry'])
        if not geom.is_valid:
            geom = geom.buffer(0)
        if geom.is_valid:
            pays[liste_pays['properties']['A3']] = geom

Ce n'est maintenant plus qu'une question de boucle :

for code, pays in liste_pays.items():
    if zone.centroid.within(pays):
        print(code)  # On a trouvé le pays !

Conclusion

Il est donc très simple de mettre en place des calculs géographique, et ce sans avoir besoin d'installer un serveur tel que PostGIS.

N'hésitez pas à consulter nos Formations SIG.