clarisma / geodesk-py

Fast and storage-efficient spatial database engine for OpenStreetMap data
https://docs.geodesk.com/python
GNU Lesser General Public License v3.0
36 stars 0 forks source link

Incorrect result for `Features.within` #57

Closed osmuser63783 closed 2 months ago

osmuser63783 commented 6 months ago

Features.within returns the following island as lying entirely inside the river area, which is wrong.

Screenshot 2024-04-28 at 17 37 12 Screenshot 2024-04-28 at 17 37 57

I think this is a bug in Geodesk.

Steps to reproduce:

my_water = next(w for w in planet_240422("a[natural=water]").ways if w.id == 854475224)
my_island = planet_240422("a[place=island,islet]").within(my_water).first # this line returns the island but shouldn't 
my_island.shape.within(my_water.shape) # Shapely returns false, which is the expected result

As a side note the definitions of within etc. can be a little counterintuitive: An island isn't "within" a lake if it shares a node with the lake, etc. If Geodesk follows the same definitions as Shapely, it could be helpful to link to this somewhere in the Geodesk docs.

clarisma commented 6 months ago

Thanks for opening this (the detailed screenshots are especially helpful). Yes, this is clearly a bug.

On the definition of within, per OGC a shape is considered to be within itself (so A can share nodes with B and still be considered within B). This is the definition used by GEOS (and hence Shapely). But in your example, the island lies outside the river area (all its nodes lie inside the river area, which appears to trigger the false positive).