On Sat, 28 Nov 2020 at 22:24, Paul Ramsey <pram...@cleverelephant.ca> wrote:
> > > On Nov 28, 2020, at 1:13 PM, Paul Ramsey <pram...@cleverelephant.ca> > wrote: > > > >> On Nov 28, 2020, at 1:11 PM, Paul Ramsey <pram...@cleverelephant.ca> > wrote: > >> > >>> On Nov 28, 2020, at 1:05 PM, Joris Van den Bossche < > jorisvandenboss...@gmail.com> wrote: > >>> > >>> On Sat, 28 Nov 2020 at 21:53, Paul Ramsey <pram...@cleverelephant.ca> > wrote: > >>> > >>>> On Nov 28, 2020, at 12:44 PM, Joris Van den Bossche < > jorisvandenboss...@gmail.com> wrote: > >>>> > >>>> Thanks for trying to reproduce it in C/C++. One obvious difference > that I can spot is that we use an integer for the "item" that gets > inserted, and not the geometry itself, but I wouldn't expect that to > influence the result. > >>>> Although, trying to update your test case to do that, the test fails. > But that might also be an issue on my side due to my limited C++ experience > (it already fails on the "geoms.size()" check): > >>> > >>> Nope still not seeing it... one small mistake in your query > >>> > >>> Also not if you try my original code *with* the mistake? > >> > >> Well, with the mistake the code is saying "here's a null pointer" > (since that's what casting 0 to a void* will get you) index this. I'll see! > It shouldn't really matter from an index point-of-view, it should still > index it and return a null pointer out the back end. > > > > And indeed it does not work. So, something to examine... why can we not > index with a null item? > > Because I was using a non-null item pointer as a cue that the node was a > leaf node, so it wasn't picking up null items as leaves. > So, I've pushed that fix, does that change your test suite all that much? > > Yes, thanks, that fixes our tests. Now, as mentioned on the commit comments on github ( https://github.com/libgeos/geos/commit/ad0122541f85886cf9eb9dd65d8246e0144f7ffa#commitcomment-44623679), if you rather not support this on the GEOS side, I am fine with having to fix this in PyGEOS to have it working correctly with GEOS 3.9. I quickly started looking at a fix in PyGEOS at https://github.com/pygeos/pygeos/pull/261 (it's not yet working though, I suppose because we need to keep track of the inserted items) > P > > > > > > P > > > >> > >>> Because your version indeed passes for me as well, but I *think* the > version I wrote resembles more closely the PyGEOS code (so it might be an > issue in our C code on how we use the tree). > >>> > >>> , trying to cast the int to a void, instead of passing in a the > pointer to the address, here's one that works. Unfortunately that leaves us > no closer to knowing why the SimpleSTRtree is unhappy in the python > context. I fear I may just have to revert the CAPI to the old tree. > >>> > >>> > >>> // querying tree with box > >>> template<> > >>> template<> > >>> void object::test<9> > >>> () > >>> { > >>> GEOSSTRtree* tree = GEOSSTRtree_create(10); > >>> > >>> GEOSGeometry* g = GEOSGeomFromWKT("POINT (2 3)"); > >>> int payload = 876; > >>> GEOSSTRtree_insert(tree, g, &payload); > >>> > >>> GEOSGeometry* q = GEOSGeomFromWKT("POLYGON ((0 0, 10 0, 10 10, 0 10, > 0 0))"); > >>> > >>> typedef std::vector<int*> IList; > >>> IList items; > >>> ensure_equals(items.size(), 0); > >>> GEOSSTRtree_query( > >>> tree, > >>> q, > >>> [](void* item, void* userdata) { > >>> IList* items = (IList*)userdata; > >>> items->push_back((int*)item); > >>> }, > >>> &items); > >>> > >>> ensure_equals(items.size(), 1); > >>> > >>> ensure_equals(*(items[0]), payload); > >>> > >>> GEOSGeom_destroy(q); > >>> GEOSGeom_destroy(g); > >>> GEOSSTRtree_destroy(tree); > >>> } > >>> > >>> > >>> > >>> > >>>> > >>>> --- a/tests/unit/capi/GEOSSTRtreeTest.cpp > >>>> +++ b/tests/unit/capi/GEOSSTRtreeTest.cpp > >>>> @@ -268,10 +268,11 @@ void object::test<8> > >>>> { > >>>> GEOSSTRtree* tree = GEOSSTRtree_create(10); > >>>> GEOSGeometry* g = GEOSGeomFromWKT("POINT (2 3)"); > >>>> - GEOSSTRtree_insert(tree, g, g); > >>>> + int idx = 0; > >>>> + GEOSSTRtree_insert(tree, g, (void*)idx); > >>>> GEOSGeometry* q = GEOSGeomFromWKT("POLYGON ((0 0, 10 0, 10 10, 0 > 10, 0 0))"); > >>>> > >>>> - typedef std::vector<GEOSGeometry*> GList; > >>>> + typedef std::vector<int> GList; > >>>> GList geoms; > >>>> ensure_equals(geoms.size(), 0); > >>>> GEOSSTRtree_query( > >>>> @@ -279,23 +280,16 @@ void object::test<8> > >>>> q, > >>>> [](void* item, void* userdata) { > >>>> GList* geoms = (GList*)userdata; > >>>> - geoms->push_back((GEOSGeometry*)item); > >>>> + geoms->push_back(*((int *)item)); > >>>> }, > >>>> &geoms); > >>>> > >>>> ensure_equals(geoms.size(), 1); > >>>> - const GEOSCoordSequence* seq = GEOSGeom_getCoordSeq(geoms[0]); > >>>> - > >>>> - double x = -1; > >>>> - double y = -1; > >>>> - GEOSCoordSeq_getXY(seq, 0, &x, &y); > >>>> - ensure_equals(x, 2.0); > >>>> - ensure_equals(y, 3.0); > >>>> + ensure_equals(geoms.at(0), 0); > >>>> > >>>> > >>>> On Sat, 28 Nov 2020 at 20:55, Paul Ramsey <pram...@cleverelephant.ca> > wrote: > >>>> > >>>> > >>>>> On Nov 28, 2020, at 8:11 AM, Joris Van den Bossche < > jorisvandenboss...@gmail.com> wrote: > >>>>> > >>>>> Hi all, > >>>>> > >>>>> On the CI of PyGEOS we have a build testing against GEOS master, and > somewhere in the last 4 days, a lot of the STRtree tests started failing > (see eg https://github.com/pygeos/pygeos/runs/1465460418#step:9:86). > Looking at the commits of the last days, this might be related to the > SimpleSTRtree work? > >>>>> > >>>>> A small (python) example of a tree consisting of a single point, > which no longer is returned when querying the tree with a big polygon that > certainly contains the point: > >>>>> > >>>>> Using released version of GEOS: > >>>>> > >>>>>>>> import pygeos > >>>>>>>> pygeos.geos_version > >>>>> (3, 8, 1) > >>>>>>>> point = pygeos.Geometry("POINT (2 3)") > >>>>>>>> tree = pygeos.STRtree([point]) > >>>>>>>> tree.query(pygeos.box(0, 0, 10, 10)) > >>>>> array([0]) > >>>>> > >>>>> This is correctly returning the index of the single point. But when > running with the latest GEOS master, the query doesn't find any point of > the tree: > >>>>> > >>>>>>>> import pygeos > >>>>>>>> pygeos.geos_version > >>>>> (3, 9, 0) > >>>>>>>> point = pygeos.Geometry("POINT (2 3)") > >>>>>>>> tree = pygeos.STRtree([point]) > >>>>>>>> tree.query(pygeos.box(0, 0, 10, 10)) > >>>>> array([], dtype=int64) > >>>>> > >>>>> Are there changes expected in how the STRtree C API functions or > required changes in user code? Or maybe we are using it in some > incorrect/unexpected way? (code is at > https://github.com/pygeos/pygeos/blob/master/src/strtree.c) > >>>> > >>>> There are changes, I don't think you're mis-using anything. I swapped > the CAPI to use the SimpleSTRtree, figuring it would be good to share the > performance win with downstream. However, I can swap it back to the > original STRtree if this remains a problem. > >>>> > >>>> One thing I noticed when trying to construct GEOS envelopes directly > was that annoyingly they were xmin xmax, ymin ymax, but I doubt that would > be a problem in your pre-existing working test. > >>>> > >>>> I just reconstructed your test in the GEOS CAPI suite, and it works > as one would expect. (Namely, it finds the one point.) So I'm not sure why > your test is getting different results. > >>>> > >>>> > >>>> // querying tree with box > >>>> template<> > >>>> template<> > >>>> void object::test<8> > >>>> () > >>>> { > >>>> GEOSSTRtree* tree = GEOSSTRtree_create(10); > >>>> GEOSGeometry* g = GEOSGeomFromWKT("POINT (2 3)"); > >>>> GEOSSTRtree_insert(tree, g, g); > >>>> GEOSGeometry* q = GEOSGeomFromWKT("POLYGON ((0 0, 10 0, 10 10, 0 > 10, 0 0))"); > >>>> > >>>> typedef std::vector<GEOSGeometry*> GList; > >>>> GList geoms; > >>>> ensure_equals(geoms.size(), 0); > >>>> GEOSSTRtree_query( > >>>> tree, > >>>> q, > >>>> [](void* item, void* userdata) { > >>>> GList* geoms = (GList*)userdata; > >>>> geoms->push_back((GEOSGeometry*)item); > >>>> }, > >>>> &geoms); > >>>> > >>>> ensure_equals(geoms.size(), 1); > >>>> const GEOSCoordSequence* seq = GEOSGeom_getCoordSeq(geoms[0]); > >>>> > >>>> double x = -1; > >>>> double y = -1; > >>>> GEOSCoordSeq_getXY(seq, 0, &x, &y); > >>>> ensure_equals(x, 2.0); > >>>> ensure_equals(y, 3.0); > >>>> > >>>> GEOSGeom_destroy(q); > >>>> GEOSGeom_destroy(g); > >>>> GEOSSTRtree_destroy(tree); > >>>> } > >>>> > >>>> > >>>> > >>>> > >>>> > >>>>> > >>>>> Best, > >>>>> Joris > >>>>> > >>>>> > >>>>> On Wed, 25 Nov 2020 at 00:44, Paul Ramsey <pram...@cleverelephant.ca> > wrote: > >>>>> Hey all, just truing up what's underway and nearly there... > >>>>> > >>>>> - Am I right that Z coordinates are nearly done? What's the status > there? > >>>>> > >>>>> I've been trying to address some performance issues, with some > success and some ... other things. > >>>>> > >>>>> The success is the SimpleSTRtree, which is just the current STRtree > but without the inheritance structure and with the nodes stored all next to > each other in contiguous memory for more locality. For at least one use > case I've seen 20% speed-ups on overlays, using the SimpleSTRtree in place > of the STRtree inside the MCIndexNoder. I have not seen any slow-downs. I > have pushed the SimpleSTRtree into master. > >>>>> > >>>>> While I have implemented the nearestNeighbor() functionality on the > SimpleSTRtree, I haven't hooked it up to anything yet. It could go into the > IndexedFacetDistance, if anyone is super enthusiastic about it. From there > it would affect searching in PreparedGeometry of various sorts. > >>>>> > >>>>> I also tried using a similar trick with the MonotoneChainBuilder > that sits inside the MCIndexNoder, replacing individual heap allocations > with slabs by putting objects onto a std::deque, and incidentally stripping > out some book-keeping. While that seems to pick up about 3-5% speedwise, > unfortunately something about my implementation is incorrect (and in a > wonderfully subtle way) as it fails testing on some platforms (not mine). > https://github.com/pramsey/geos/tree/monotone-chain-builder > >>>>> > >>>>> I've put that work to the side for now. > >>>>> > >>>>> All the performance talk is mostly because JTS still runs a lot > faster than GEOS for some bulk processing. My current test is a big union > of watershed boundaries, about 6MB of data, which takes about 20s under > GEOS and about 25% of that under JTS. It's a big gap, and in theory the two > code bases are pretty aligned right now. Same overlayNG engine, etc. So I > figure there has to be a big implementation ball of performance hiding > under the covers somewhere. No luck thus far. > >>>>> > >>>>> I think we're close, looking forward to release :) > >>>>> > >>>>> P > >>>>> _______________________________________________ > >>>>> geos-devel mailing list > >>>>> geos-devel@lists.osgeo.org > >>>>> https://lists.osgeo.org/mailman/listinfo/geos-devel > >>>>> _______________________________________________ > >>>>> geos-devel mailing list > >>>>> geos-devel@lists.osgeo.org > >>>>> https://lists.osgeo.org/mailman/listinfo/geos-devel > >>>> > >>>> _______________________________________________ > >>>> geos-devel mailing list > >>>> geos-devel@lists.osgeo.org > >>>> https://lists.osgeo.org/mailman/listinfo/geos-devel > >>>> _______________________________________________ > >>>> geos-devel mailing list > >>>> geos-devel@lists.osgeo.org > >>>> https://lists.osgeo.org/mailman/listinfo/geos-devel > >>> > >>> _______________________________________________ > >>> geos-devel mailing list > >>> geos-devel@lists.osgeo.org > >>> https://lists.osgeo.org/mailman/listinfo/geos-devel > >>> _______________________________________________ > >>> geos-devel mailing list > >>> geos-devel@lists.osgeo.org > >>> https://lists.osgeo.org/mailman/listinfo/geos-devel > > _______________________________________________ > geos-devel mailing list > geos-devel@lists.osgeo.org > https://lists.osgeo.org/mailman/listinfo/geos-devel >
_______________________________________________ geos-devel mailing list geos-devel@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/geos-devel