4

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.

PGAdmin Geometry Viewer

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:

Map on epsg.io

They contradict themselves!

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

5
  • 4
    It is either lat-lon or lon-lat and it depends on the software. Sometimes you can even select it with an option. Also the format may have an effect, see GML and GML 3 examples in postgis.net/docs/ST_AsGML.html. For PostGIS WKT is supposed to be in lon-lat order. Somebody will say that lat-lon is right and lon-lat is wrong but it does not help. Both are used and you just need to know when to use which order.
    – user30184
    Jun 20 at 10:31
  • this is horrible 😮
    – Klamann
    Jun 20 at 11:02
  • What you are inserting is a WKT description of a point, and if you look into the description of WKT point, it states: "Coordinates for geometries may be 2D (x, y)" - no matter which CRS is used, so first east - west, then north south. It may seem a bit unusual for latlon, but having it switch back and forth if it was latlon or a projected system would have been more difficult to handle. Jun 20 at 11:03
  • 3
    Unfortunately it is not defined that x means longitude/easting and y means latitude/northing. They are just x and y but very many who implemented WKT did it that way. Some have changed their minds later like dev.mysql.com/doc/relnotes/mysql/8.0/en/news-8-0-1.html Incompatible Change: These functions now interpret latitude and longitude coordinates as in the order specified by the spatial reference system. The functions also accept an optional argument to override the default axis order.
    – user30184
    Jun 20 at 11:09
  • I think that SRID=4326 isn't EPSG:4326 but OCG:CRS84. Treat it that way. Jun 22 at 11:53

2 Answers 2

4

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:

6
  • thanks for the clarification! So it seems like PostGIS is "wrong" here, or at least not quite consistent with EPSG:4326. Is there a good way to deal with this in PostGIS? Some option that allows me to specify latitude first?
    – Klamann
    Jun 20 at 11:26
  • 2
  • I'm not sure that it always has - the WMS 1.0 spec certainly defined it as lon,lat (which is one of the problems)
    – Ian Turton
    Jun 20 at 14:43
  • about ST_FlipCoordinates, when would I apply this? Should I store the data lat,lon in the database and then flip whenever I want to visualize it? Or store it lon,lat and flip whenever I want to extract it? What if someone forgets to flip the coordinates when they query the data? I don't know, this seems like a setup for failure
    – Klamann
    Jun 20 at 15:30
  • 2
    Agree that OGC confused matters by incorrectly swapping the axis order in early standards, see directionsmag.com/article/2398
    – nmtoken
    Jun 20 at 16:14
1

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.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.