Long story short, you create a table with a geometry column that stores points with srid 4326 (read geographic lat/lng, or should I write lng/lat, coordinates). That should be interpreted as geography right? Not at all!
To demonstrate, we’ll need to use some GIS related functions, as using pg_typeof only yields “unknown.
db=> SELECT pg_typeof('POINT(45.5 -73.5)');
pg_typeof
-----------
unknown
Let’s use the ST_Distance function and have a look at its documentation.
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. Ifuse_spheroidis false, a faster spherical calculation is used.
Ok, so there are two function signatures. Let’s start with the first one taking two geometries as input.
db=> SELECT ST_Distance('POINT(45.5 -73.5)', 'POINT(45 -73.5)');
st_distance
-------------
0.5
Two geometries, the distance is calculated in degres: 0.5, all good, nothing worth discussing here.
Let’s try the second function signature by setting the points as srid 4326, and call the function again.
db=> SELECT ST_Distance('SRID=4326;POINT(45.5 -73.5)', 'SRID=4326;POINT(45 -73.5)');
st_distance
-------------
0.5
Wait what? Not the result one may expect. So, as the title of this post implies, setting SRID to 4326 is not enough, it’s still a geometry. To have a proper geography, the values need to be casted.
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! The answer we were looking for. Ok, I’ve learned something. If you want to keep going, here is another “fun” thing I’ve discovered. While you cannot mix SRID within a query, you can apparently mix geometries and geographies. Two more examples. The firsts one, mixing SRIDs does not work.
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
The second one, mixing geometries and geographies works, converting any of the non geographies one to geographies.
db=> SELECT ST_Distance('SRID=4326;POINT(45.5 -73.5)', CAST('SRID=4326;POINT(45 -73.5)' as geography));
st_distance
---------------
15857.0461096
That’s not very consistent.
Finally, for any one trying to convert geometries to geographies in sqlalchemy, you need to import the cast function, you need the geoalchemy2 library to import its Geography type. Then you cast the field to a Geography by setting the arguments to None. Otherwise it produces a default conversion with SRID=-1, and that does not work. Here is an example over a field.
from geoalchemy2.types import Geography
from sqlalchemy.sql import cast
# assuming a model.field exists
field_as_geography = cast(model.field, Geography(None))