duckdb / duckdb_spatial

MIT License
483 stars 37 forks source link

[BUG][0.10.0] ST_Distance_Spheroid() seems to be wrong #287

Closed keltia closed 7 months ago

keltia commented 7 months ago

I have a Rust program to calculate geodesic distances with both my implementations and the geo crate. Distances are pretty similar but as soon as I use ST_Distance_Spheroid(), I get wildly different results.

Code is at: https://github.com/keltia/fetiche-rs/blob/develop/process-data/examples/dist.rs

Results are:

warning: `process-data` (example "dist") generated 2 warnings (run `cargo fix --example "dist"` to apply 2 suggestions)
    Finished dev [unoptimized + debuginfo] target(s) in 5.36s
     Running `/Users/roberto/Src/Rust/src/fetiche-rs/target/debug/examples/dist 49.616426 6.174652 49.6300506592 6.2211528048`
Distance between
  Point(Coord { x: 6.174652, y: 49.616426 })
and
  Point(Coord { x: 6.2211528048, y: 49.6300506592 })
Distances:
3676.33 m haversines
3676.33 m (sin/cos)
3685.81 m geo::geodesic
3676.29 m geo::haversines
3685.81 m geo::vincenty
5358.91 m duckdb:speh

Did I miss something? The functions are still undocumented so I don't know.

The query is basic:

    let dbh = duckdb::Connection::open_in_memory()?;
    dbh.execute("LOAD spatial", [])?;

    let distduck: f64 = dbh.query_row("SELECT ST_Distance_Spheroid(ST_Point(?, ?), ST_Point(?, ?))",
                                      params![point1.longitude, point1.latitude, point2.longitude, point2.latitude], |row| {
            Ok(row.get_unwrap(0))
        },
    )?;
Maxxen commented 7 months ago

Hi! Thanks for reporting this issue!

The _spheroid functions expect WGS84/EPSG:4326 which mandates a latitude/longitude axis order, not a longitude/latitude axis order, which is why you are getting wrong result. Admittedly this is somewhat confusing and made worse by the fact that DuckDB spatial is in the minority of trying to remain ambivalent and respect the axis order mandated by the projection compared to most other GIS tools that always do lon/lat. This will probably change on our end in the future.

Maxxen commented 7 months ago

I have some work in progress on new docs that explain this here

keltia commented 7 months ago

Thanks for the answer. This also points to a bigger issue, there is no documentation whatsoever on these functions... 😅 I sent a comment on the now-closed PR #74 before 0.10.0.

keltia commented 7 months ago

Much better, thanks!

     Running `C:\Users\roberto\Src\rs\src\fetiche-rs\target\debug\examples\dist.exe 49.616426 6.174652 49.6300506592 6.2211528048`
Distance between
  Point(Coord { x: 6.174652, y: 49.616426 })
and
  Point(Coord { x: 6.2211528048, y: 49.6300506592 })
Distances:
3676.33 m haversines
3676.33 m (sin/cos)
3685.81 m geo::geodesic
3676.29 m geo::haversines
3685.81 m geo::vincenty
3685.81 m duckdb:speh
keltia commented 7 months ago

I have some work in progress on new docs that explain this here

Ok, I see it now. Thanks.