Feature: Spherical calculations

[Maxxen]

While I'm hesitant to implement a full blown GEOGRAPHY type, I do think that having a couple of spatial functions using geodetic algorithms operating on a sphere would be a good complexity/usefulness compromise.

• ST_Geodetic_Area
• ST_Geodetic_Length
• ST_Geodetic_Distance

[blackrez]

I don't if it is help, but actually I used the rust port of https://geographiclib.sourceforge.io/C++/doc/index.html for geodesic calculation.

[bfgdigital]

I'm currently trying to replicate the following

ST_DistanceSphere(
        ST_MakePoint(lon, lat),
        ST_MakePoint(last_lon, last_lat)
 ) AS distance_m,

with

ST_Distance(
        ST_Point(lon, lat),
        ST_Point(last_lon, last_lat)
) AS distance_m,

but not being a spatial guy, I'm unsure of the implications of this.
From what I've discerned I would require these Spherical calculations to provide the 6370986 meter radius effect on the measurements between the two points.

this sounds like a pretty fundamental requirement of geo data so I'm adding this here in hope that it isn't required, but perhaps does require spherical calculations.

[Maxxen]

To clarify, all the algorithms implemented as of now assume a planar coordinate system. This makes a lot of algorithms much easier and faster to implement and works well when you are processing geospatial data that lies within a small-ish region, say e.g. on the scale of a country or two. Once you reach a scale where the curvature of the earth becomes important, then yeah, measurements like distances and areas become extremely inexact if implemented with straight forward cartesian geometry . The PostGIS documentation has a great example of this.

Note that in general we don't really know the exact shape of the earth, so there is really three different ways to approximate geospatial calculations with progressively increasing difficulty, performance overhead and precision.

1. Project into a planar coordinate space. This is what we currently achieve with GEOS and PROJ.
2. Model the earth as a perfect sphere. This is what e.g. Googles S2 library does (and I think BigQuery?).
3. Model the earth as an ellipsoid. This is what the algorithms in GeographicLib does.

Like PostGIS it would be cool if we could support all three to some degree, but a lot of the more advanced algorithms (e.g. polygon intersection) are already extremely complicated in the planar case, so I'd like to take this in small steps.

That said, I have something in the works using Geographiclib for the elementary algorithms (distance, length, area). But I think we should implement and focus on 2. in the future for a potential GEOGRAPHY type.

[Maxxen]

Partially implemented in #74

[keltia]

Guys guys... #74 was merged in May 2023 but the documentation was never updated so I just discovered that they are in the current extension!!!

D select ST_Distance_spheroid(ST_Point(49.60, 6.12), ST_Point(49.60, 6.14));
┌────────────────────────────────────────────────────────────────────┐
│ st_distance_spheroid(st_point(49.60, 6.12), st_point(49.60, 6.14)) │
│                               double                               │
├────────────────────────────────────────────────────────────────────┤
│                                                  1445.776772642289 │
└────────────────────────────────────────────────────────────────────┘

[sonjoonho]

Hi – it's exciting to see that this extension is still being actively developed! Almost a year on since this roadmap item was opened, I was wondering if the calculation of geographic/spheroid distances between polygons and points was on the horizon? Otherwise it looks like falling back to PostGIS is my best bet (but open to other suggestions).

Line segment <-> point and the additional combinations that follow require an iterative algorithm that I've left out for now.

I would also be curious to know what this algorithm involves. Thanks!

[Maxxen]

Time flies!
Ive recently completed a pretty major overhaul of the internal geometry representation used during execution so im planning on going back and adding the remaining missing measurement functions, including geodetic/haversine based ops, in the near future.

Im currently working on refactoring the docs and will try to shape up and share more about the roadmap for my plans for 2024 as well.

[jtbaker]

DuckDB version == v1.0.0

Thanks for your contributions on this! I think I've identified a possible issue where the x and y coordinates are getting transposed for each other in one of the checks. If you go past -90 (min latitude but not min longitude) on the x coordinate (longitudinal) plane, the calculations seem to stop working

SELECT ST_AREA_SPHEROID(ST_ENVELOPE(ST_BUFFER(ST_POINT(-89.999, 30), .0001)));

st_area_spheroid(st_envelope(st_buffer(st_point(-89.999, 30), .0001))) │
│                                 double                                 │
├────────────────────────────────────────────────────────────────────────┤
│                                                   0.008709549903869629 

^ this works as expected, while incrementing up to from -89.999 to -90 for the longitudinal coordinate returns a nan with everything else the same

SELECT ST_AREA_SPHEROID(ST_ENVELOPE(ST_BUFFER(ST_POINT(-90, 30), .0001)));

 st_area_spheroid(st_envelope(st_buffer(st_point(-90, 30), .0001))) │
│                               double                               │
├────────────────────────────────────────────────────────────────────┤
│                                                                nan 

[mtravis]

@Maxxen I'm also seeing the error above.

overturemaps download -f geoparquet --bbox -71.063432,42.301047,-71.009501,42.382221 -o boston.parquet -t building

D SELECT
  id,
  ST_AREA(ST_GeomFromWKB(geometry)) as area,
  ST_Area_Spheroid(ST_geomFromWKB(geometry)) AS cartesian_area
  FROM 'boston.parquet' limit 5;
┌──────────────────────────────────┬────────────────────────┬───────────────────┐
│                id                │          area          │  cartesian_area   │
│             varchar              │         double         │      double       │
├──────────────────────────────────┼────────────────────────┼───────────────────┤
│ 08b2a3066bd41fff0200fc6ab276c90d │  1.285837500012519e-08 │  51.9845110103488 │
│ 08b2a3066bd4cfff0200bc057338d4ed │ 1.2229789998652358e-08 │ 49.44350370578468 │
│ 08b2a3066bd41fff020028311e9c371c │ 1.5778415000042524e-08 │  63.7901326790452 │
│ 08b2a3066bd4cfff02008dcf88943bb9 │ 1.1989005000108085e-08 │ 48.47045111026091 │
│ 08b2a3066bd4cfff020002be80094552 │ 2.0510755000441452e-08 │ 82.92415748536587 │

overturemaps download -f geoparquet --bbox -122.365604,47.573155,-122.307926,47.623162 -o seattle.parquet -t building

D SELECT
  id,
  ST_AREA(ST_GeomFromWKB(geometry)) as area,
  ST_Area_Spheroid(ST_geomFromWKB(geometry)) AS cartesian_area
  FROM 'seattle.parquet' limit 5;
┌──────────────────────────────────┬────────────────────────┬────────────────┐
│                id                │          area          │ cartesian_area │
│             varchar              │         double         │     double     │
├──────────────────────────────────┼────────────────────────┼────────────────┤
│ 08b28d555096bfff0200ba6e244e17a4 │   8.85314630001753e-07 │            nan │
│ 08b28d5550948fff0200b8ce7ac64e3b │ 1.0191719999341135e-08 │            nan │
│ 08b28d555095efff02004bd14aeeb6be │  9.111243699996493e-07 │            nan │
│ 08b28d5550958fff0200eb9c47a75fc9 │  2.090709999588343e-09 │            nan │
│ 08b28d5550855fff02007c80ff8ad9f0 │  2.003438000002535e-08 │            nan │
└──────────────────────────────────┴────────────────────────┴────────────────┘

[kentstephen]

I just want to weigh in here that I believe this is the root of an issue I am having with Overture Maps, specifically with:

ST_Transform(ST_GeomFromWKB(geometry), 'EPSG:4326', 'EPSG:3857') as reprojected_geometry

I then calculate the ST_Area(reprojected_geometry)

When I run that on the whole world I get NULL values for a good chunk of it:

I am a fairly novice, self taught developer - I don't understand the specific details of this kind of thing, but it's my understanding that ST_Transform uses spherical calculations.

[Maxxen]

Hi! There's a couple of things at play here.

ST_Transform defaults to respecting the axis order defined by the CRS authority. EPSG:4326 is defined as lat/lon, while 'EPSG:3857' is defined as lon/lat. Meaning that when transforming 'EPSG:4326 to EPSG:3857 it is expected that the input is already in lat/lon axis order.

You can pass a boolean as an optional third argument to force st_transform to always expect [lon/lat] axis order regardless of what projection you use. You can also use st_flipcoordinates to adjust the axis order manually. The new docs are still a bit rough around the edges but I've attempted to document this better here: . https://github.com/duckdb/duckdb_spatial/blob/main/docs/functions.md#st_transform

Additionally, while most functions in the spatial extension is ambivalent to the axis order or projection (as they almost all work in planar space) the _spheroid functions expect their input to be in lat/lon EPSG:4326 (WGS84).

This is a common cause of confusion, see e.g:
#222, #219 #215 #179 #63

I've elaborated more on this both on discord and in the issues reported above. But in short:

Most GIS tools tend to always work with a fixed axis order. Most FOSS GIS tools prefer lon/lat (x/y) while others (especially in the GPS/radio-space) prefer lat/lon. Around the time DuckDB spatial was created it seemed like to me that the big libraries in the open source GIS world such as GDAL and PROJ had recently made an attempt to actually start respecting the CRS authority instead of always expecting lon/lat. See e.g. GDAL axis order post 3.0 and PROJ's FAQ. Therefore the approach taken in DuckDB spatial was to simply stay ambivalent to the axis order used. As we don't store SRID's in each geometry anyway it made sense to me that the user would ultimately be responsible for knowing what coordinate system and axis order they worked within.

Unfortunately it seems like the idea of respecting the axis order defined by the CRS authority never really caught on in practice. Modern specs such as GeoJSON uses geographic coordinates but mandates lon/lat anyway, and even brand new specs such as GeoParquet require lon/lat always.

Judging by the number of issues opened by users confused by this, I am definitely ready to change this behavior so that we also normalize all coordinates into easting/northing (e.g. lon/lat) when ingested by duckdb internally. However, this will require a bit of work as we then would need some way to track the projection information for each geometry column whenever we import/export data from duckdb spatial so that we can flip the coordinates back from/to what the CRS mandates when data enters or leaves spatial.

I have some ideas on how to do this but it is going to be part of a larger overhaul I've got planned for how we handle projection information in spatial, so I don't have a timeline for when this behavior will change.