duckdb_spatial

duckdb_spatial copied to clipboard

When reading geometries won't be loaded properly and calculations will be wrong if the order of lat/long isn't correct. So flip_coordinates should be used. In the short term, maybe give it its own section in the readme/documentation? (right now it's mentioned in passing for the GeoJSONSeq write). Longer term, the {sf} package in R handles it perfectly. I don't know how.

More generally, PostGIS is long/lat for constructors (https://postgis.net/docs/ST_Point.html). Apache Sedona had the same issue and ended up changing to match PostGIS (https://github.com/apache/sedona/pull/639).

For example this dataset of railways: https://geodata.bts.gov/datasets/d83e85154a304da995837889cc4012e3_0/about 'KM' has the right length.

CREATE TABLE rail as SELECT *, st_geomfromwkb(wkb_geometry) as geometry
                       FROM st_read('./North_American_Rail_Network_Lines/North_American_Rail_Network_Lines.shp')
                       WHERE tracks > 0


 SELECT
    KM,
    st_length(st_transform(st_flipcoordinates(geometry), 'EPSG:4326', 'EPSG:5070')) AS test_len_right,
    st_length(st_transform(geometry, 'EPSG:4326', 'EPSG:5070')) AS test_len_wrong
  FROM rail_2
  WHERE (STATEAB = 'NY')

Hi! The GEOMETRY is internally ambivalent to which axis order is used, as far as we know there's only an "x" and a "y", coordinate, whatever corresponds to lat or long is completely dependent on what data you feed the extension, or what CRS you work with.

I realise however that this may be problematic once we introduce geodesic functions that do assume a specific layout. In the case of reading data properly I don't know what capabilities GDAL has to extract this information, but I suspect it may be possible to get it from PROJ if the CRS is known at import. I'll look into how {sf} handles this in the future.

Stumbled upon this too. It usually is the case with WGS-84/EPSG:4326

Currently very confusing:

PostGIS

SELECT
  ST_GeomFromText('POINT(13 52)', 4326) -- SRID=4326;POINT(13 52)
  , ST_Transform(ST_GeomFromText('POINT(13 52)', 4326), 25833) -- SRID=25833;POINT(362705.634104228 5762926.812790221)
;

DuckDB

SELECT 
  ST_GeomFromText('POINT(13 52)') -- POINT (13 52)
  , ST_Transform(ST_GeomFromText('POINT(13 52)'),'EPSG:4326','EPSG:25833') -- POINT (4788172.580893706 1785584.7956043417)
  , ST_Transform(ST_FlipCoordinates(ST_GeomFromText('POINT(13 52)')),'EPSG:4326','EPSG:25833') -- POINT (362705.634104228 5762926.812790221)
;

I hear you. I'll work on aligning with what PostGIS does (always lon/lat, ignore CRS axis-order definition) next chance I get.

Thanks @Maxxen. Great work by the way, sorry to not mention this before. 😬

Hi, I see this problem is already open, the current behavior is really confusing, if you have the ability to match it with how the transformation footer behaves with other databases like Postgis I would be very grateful, it's a great extension.

You might want to apply this to PJ objects... https://proj.org/en/9.4/development/reference/functions.html#c.proj_normalize_for_visualization

Hello, and thank you for this awesome extension to the awesome DuckDB! 😄

I'm coming from the PostGIS world and wanted to drop in a +1 that this particular behavior was confusing to me. I was I was doing a spatial join with an ST_Transform() in the middle and getting no results. After some debugging I realized the Xs and Ys were flipped which led me to this issue. I've never had to use PostGIS ST_FlipCoordinates() — in fact, I didn't know it existed. PostGIS must automatically flip coordinates when appropriate.

If there's a project goal of reducing friction for others coming from tools like PostGIS, I would strongly suggest reconsidering this as a default — just my two cents!