Order of latitude and longitude in EPSG:4326

Question

I recently started working with GPS data in PostGIS and am very confused about the order of latitude and longitude. I'm from Europe, and here everyone I know uses latitude first, longitude second, e.g. (48.137, 11.575) points to Munich (and not somewhere off the coast of Somalia). This is also what Google Maps, OpenStreetMap and many other services will tell me, when I enter these coordinates in the search box.

Now I created a table in PostGIS

CREATE TABLE public.epsg_4326 (
    id integer primary key,
    geom geometry(Point, 4326)
);

The manual I followed said EPSG:4326 is basically WGS84, which is used by GPS, which is what I want, so I picked it for my geometry type.

Now let's insert some data...

INSERT INTO public.epsg_4326 (id, geom) VALUES (1, 'POINT(48.137 11.575)');

Then I inspected the result in PGAdmin. When I switch to the geometry viewer, it will place my point in Africa. So clearly, PostGIS (or pgAdmin) expect my data as long/lat, instead of lat/long.

I tried to figure out how EPSG:4326 is defined and at https://epsg.io/4326 they clearly state

Coordinate system: Ellipsoidal 2D CS. Axes: latitude, longitude. Orientations: north, east.

But, if I follow the link to the map on the very same page, this is what happens:

They contradict themselves!

What's going on here? How is EPSG:4326 really defined? Should I use latitude or longitude first?

Answer 1

Note EPSG.io despite its name is not the official registry of EPSG. The official EPSG Geodetic Parameter Dataset is https://epsg.org/home.html.

The CRS EPSG:4326 is (and always has been) defined as first axis latitude second axis longitude.

Most EPSG geographic coordinate systems follow the latitude then longitude, coordinate axis order, following the default of ISO 6709 Standard representation of geographic point location by coordinates.

Also see:

• Why latitude/longitude and not the other way round?
• Displaying coordinates and inputs as LatLon or LonLat?

Answer 2

I couldn't quite accept that PostGIS forces me to use (lon,lat) order when EPSG:4326 says it should be (lat,lon), so I spent way too much time thinking of a clever solution to find a suitable abstraction and this is what I came up with:

• On the client side, I'm working with Python, using SQLAlchemy and GeoAlchemy.
• I specified a custom GPS type that is an extension of Geometry which automatically swaps lat,lon to lon,lat on insert, and lon,lat back to lat,lon on select.
• Additionally, it can parse WKT data to Shapely shapes and limits the precision of coordinates with ST_QuantizeCoordinates to a reasonable amount.

Here's the GPS type definition

from geoalchemy2 import Geometry, functions as geo_func
from geoalchemy2.shape import to_shape
from sqlalchemy import TypeDecorator, func


class GPS(TypeDecorator):
    """
    Extension of PostGIS type 'Geometry' with the following additions:

    * Flips the order of latitude and longitude, so that from the outside, 
      we get lat,lon ordering, as it should be, but in the database, 
      it is lon,lat so that PostGIS is happy
    * limits the precision of all stored coordinates to a reasonable amount
      (no raw floats needed, reduces table size when TOAST compression 
      is available)

    :param srid: a custom SRID (defaults to EPSG:4326 for GPS)
    :param precision: number of significant decimal digits after 
           the decimal point
    :param kwargs: other parameters for the Geometry type
    """
    impl = Geometry
    cache_ok = True

    def __init__(self, srid=4326, precision=7, **kwargs):
        super().__init__(srid=srid, **kwargs)
        self.precision = precision

    def column_expression(self, col):
        # this function will be executed by the database on select
        # flip coordinates from lon,lat back to lat,lon
        return getattr(func, self.impl.as_binary)(
            geo_func.ST_FlipCoordinates(col), type_=self)

    def bind_expression(self, bindvalue):
        # this function will be executed by the database on insert
        # flip coordinates from lat,lon to lon,lat and limit precision
        return geo_func.ST_QuantizeCoordinates(
            geo_func.ST_FlipCoordinates(self.impl.bind_expression(bindvalue)),
            self.precision, type_=self)

    def process_result_value(self, value, dialect):
        # this function will be executed locally on select
        # parse the raw WKB value to a shapely geometry
        return to_shape(value)

And this is how you can use

from sqlalchemy import Column, Integer
from sqlalchemy.orm import declarative_base

Base = declarative_base()


class GeoTest(Base):
    __tablename__ = 'epsg_4326'
    id = Column(Integer, primary=True)
    position = Column(GPS(geometry_type='POINT', precision=5))

Let's try it

from shapely.geometry import Point
from sqlalchemy import create_engine
from sqlalchemy.orm import sessionmaker

engine = create_engine('postgresql://user:secret@localhost:5432/test_db', echo=True)
Base.metadata.create_all(engine)
session = sessionmaker(bind=engine)

gps = GeoTest(id=0, position=Point(48.137, 11.575).wkt)
session.add(gps)
session.commit()
# stored as wkt "POINT(11.575 48.137)" in the database
gps2 = session.query(GeoTest).first()
# restored as shapely Point(48.137, 11.575) in python
# exact float values can be different, since we are rounding in binary

I'm sure there are some edge cases that I've missed, but I'm ok with this solution for now.