C’est un piège! Postgis Geometry avec un SRID 4326 n’est pas un type Geography

Longue histoire courte, vous créez une table avec une colonne geometry qui stocke des points avec srid 4326 (lire coordonnées géographiques lat/lng, ou devrais-je écrire lng/lat). Cela devrait être interprété comme geography, n’est-ce pas? Pas du tout!

Pour le démontrer, nous devrons utiliser des fonctions liées au SIG (GIS en anglais), car l’utilisation de pg_typeof n’aboutit qu’à « unknown ».

db=> SELECT pg_typeof('POINT(45.5 -73.5)');
 pg_typeof
-----------
 unknown

Utilisons la fonction ST_Distance et jetons un coup d’œil à sa documentation (en anglais).

float ST_Distance(geometry g1, geometry g2);
For geometry types returns the minimum 2D Cartesian (planar) distance between two geometries, in projected units (spatial ref units).

float ST_Distance(geography geog1, geography geog2, boolean use_spheroid = true);
For geography types defaults to return the minimum geodesic distance between two geographies in meters, compute on the spheroid determined by the SRID. If use_spheroid is false, a faster spherical calculation is used.

Ok, il y a donc deux signatures de fonctions. Commençons par la première qui prend deux géométries en paramètres.

db=> SELECT ST_Distance('POINT(45.5 -73.5)', 'POINT(45 -73.5)');
 st_distance
-------------
         0.5

Deux géométries, la distance est calculée en degrés : 0.5, tout va bien, il n’y a rien à signaler ici.

Essayons la deuxième signature de la fonction en définissant les points comme srid 4326, et appelons à nouveau la fonction.

db=> SELECT ST_Distance('SRID=4326;POINT(45.5 -73.5)', 'SRID=4326;POINT(45 -73.5)');
 st_distance
-------------
         0.5

Pardon quoi? Ce n’est pas le résultat auquel on peut s’attendre. Ainsi, comme le titre de ce billet l’indique, il ne suffit pas d’attribuer la valeur 4326 au SRID, il s’agit toujours d’une géométrie. Pour avoir une géographie correcte, les valeurs doivent être transtypés (castées en bon québécois).

db=> SELECT ST_Distance(CAST('SRID=4326;POINT(45.5 -73.5)' AS geography), CAST('SRID=4326;POINT(45 -73.5)' AS geography));
  st_distance
---------------
 15857.0461096

Tada ! ~15,857 km ! La réponse que nous cherchions. Ok, j’ai appris quelque chose. Si vous voulez continuer, voici une autre chose « amusante » que j’ai découverte. Alors que vous ne pouvez pas mélanger les SRID dans une requête, vous pouvez apparemment mélanger les géométries et les géographies. Deux autres exemples. Dans le premier cas, le mélange des SRID ne fonctionne pas.

db=> SELECT ST_Distance('SRID=4326;POINT(45.5 -73.5)', 'POINT(45 -73.5)');
ERROR:  ST_Distance: Operation on mixed SRID geometries (Point, 4326) != (Point, 0)
CONTEXTE : SQL function "st_distance" statement 1

Le second, qui consiste à mélanger les géométries et les géographies, fonctionne en convertissant les éléments non géographiques en géographies.

db=> SELECT ST_Distance('SRID=4326;POINT(45.5 -73.5)', CAST('SRID=4326;POINT(45 -73.5)' as geography));
  st_distance
---------------
 15857.0461096

Ce n’est pas très cohérent.

Enfin, pour tous ceux qui essaient de convertir des géométries en géographies dans sqlalchemy, vous devez importer la fonction cast et vous avez besoin de la bibliothèque geoalchemy2 pour importer son type Geography. Ensuite, il faut convertir le champ en géographie en définissant l’argument sur None. Sinon, la fonction produit une conversion par défaut avec SRID=-1, et cela ne fonctionne pas. Voici un exemple minimaliste.

from geoalchemy2.types import Geography
from sqlalchemy.sql import cast

# assuming a model.field exists
field_as_geography = cast(model.field, Geography(None))

Références

Leave a Reply

Ce site utilise Akismet pour réduire le pourriel. En savoir plus sur comment les données de vos commentaires sont utilisées.