duckdb / duckdb_spatial

MIT License
467 stars 34 forks source link

ST_DWithin_Spheroid and ST_Distance_Spheroid documentation #315

Closed nshiab closed 4 days ago

nshiab commented 4 months ago

Hi!

For ST_Distance_Spheroid, the documentation says this:

Returns the distance between two geometries in meters using a ellipsoidal model of the earths surface. The input geometry is assumed to be in the EPSG:4326 coordinate system (WGS84), with [latitude, longitude] axis order and the distance is returned in meters.

But for ST_DWithin_Spheroid, the documentation doesn't specify if it's the same thing.

Is it? Should we only use ST_DWithin_Spheroid with EPSG:4326 and [latitude, longitude] axis order with the target distance in meters?

Thanks!

nshiab commented 4 months ago

Here's an example explaining my confusion.

In this query, I fetch two geojson files. The first one in the table mtl is just a random point in Montreal. The second one is three points in Montreal, Toronto, and Vancouver. I flip the coordinates and then compute the spheroid distance between the random point in Montreal and the three other points.

LOAD spatial;
LOAD https;
CREATE TABLE mtl as SELECT geom as geomMtl FROM ST_Read('https://raw.githubusercontent.com/nshiab/simple-data-analysis/main/test/geodata/files/point.json');
UPDATE mtl SET "geomMtl" = ST_FlipCoordinates("geomMtl");
CREATE TABLE cities as SELECT name, geom as geomCities FROM ST_Read('https://raw.githubusercontent.com/nshiab/simple-data-analysis/main/test/geodata/files/coordinates.geojson');
UPDATE cities SET "geomCities" = ST_FlipCoordinates("geomCities");
CREATE OR REPLACE TABLE cities AS SELECT cities.*, mtl.* FROM cities CROSS JOIN mtl;
ALTER TABLE cities ADD "distToMtl" DOUBLE; UPDATE cities SET "distToMtl" = ST_Distance_Spheroid("geomMtl", "geomCities");
SELECT * from cities;

The distances look correct.

Screenshot 2024-05-09 at 9 11 44 AM

But if I use ST_DWithin_Spheroidto filter points within 500 km of the random point in Montreal, the results are weird. I would expect to see toronto and montreal, but not vancouver!

SELECT * FROM cities WHERE ST_DWithin_Spheroid("geomCities", "geomMtl", 500000);
Screenshot 2024-05-09 at 9 14 55 AM

Any idea what's going wrong?

Thank you. 🙏

CGenie commented 1 week ago

I'm confused as well. I think ST_DWithin_Spheroid should return a boolean and not a double similar to ST_DWithin and to what PostGIS does: https://postgis.net/docs/ST_DWithin.html?

Maxxen commented 1 week ago

@CGenie you're right, that's an issue indeed. The internal function signature is wrong even though I think the function is implemented correctly. I'll push a fix tomorrow.

nshiab commented 3 days ago

Thanks, @Maxxen!