questions about RECTANGLE_INTERSECTION flag

[Martin Davis]

Excellent.

It seems there might be some bugs lurking in GEOSClipByRect.  If you come
across any failure cases or things that don't work as expected, please file
an issue in GEOS or PostGIS for them.  It would be nice to get rectangle
clipping working error-free and as fast as possible.

On Wed, Nov 27, 2024 at 9:32 AM arno <arno at renevier.net> wrote:

> Thanks a lot. I was not aware of ST_ClipByBox2D, but it seems to be what I
> need.
>
> On Tue, Nov 26, 2024 at 07:56:26AM -0800, Martin Davis wrote:
> > One reason RectangleIntersection is not enabled in the Intersection code
> > path is that it has different semantics than the OverlayNG algorithm.
> This
> > might be undesirable or confusing in some use cases.
> >
> > It is exposed by GEOSClipByRect.  And this is used by PostGIS (at least
> in
> > some circumstances).  So have you tried using ST_ClipByBox2D?
> >
> > On Tue, Nov 26, 2024 at 1:02 AM arno <arno at renevier.net> wrote:
> >
> > > Hi,
> > >
> > > Long story short:
> > > Enabling USE_RECTANGLE_INTERSECTION flag speeds up postgis
> ST_Intersection
> > > by 60x for my use case
> > >
> > > Long story long:
> > > I recently wrote a strava like heatmap. It requires computing the
> > > intersections of my gpx tracks with tile layers. The PostGIS query
> used is:
> > >
> > > ```
> > > const query = `SELECT ST_ASGeoJson(ST_Intersection(geom,
> > > ST_MakeEnvelope($1, $2, $3, $4, 4326))) AS geom from
> ${this.#tableName}`
> > > ```
> > >
> > > It is the performance bottleneck on my application. I have about 2000
> gpx
> > > tracks. On a large tile, where 1300 geometries intersect, the query
> takes
> > > about 2500ms to run. I tried various query optimizations without
> sucess.
> > >
> > > But if I rebuild libgeos with USE_RECTANGLE_INTERSECTION flag set, the
> > > query succeeds in about 40ms. (ie: 60x speedup!) (
> > > https://github.com/libgeos/geos/blob/main/src/geom/Geometry.cpp#L76)
> > >
> > > Is there a way this optimization could be enabled in libgeos?
> > >
> > > I understand there are differences between the current intersection
> > > operation and the rectangle intersection optimization:
> > >
> > > + If the geometries don't intersect, the result is an empty geometry of
> > > type geometryA. With RectangleIntersection, it is an empty
> > > GEOMETRYCOLLECTION
> > > ```
> > > % geosop -a 'LINESTRING(10 10, 20 20)' -b 'POLYGON((110 10, 120 10, 120
> > > 20, 110 20, 110 10))' clipRect
> > > GEOMETRYCOLLECTION EMPTY
> > > % geosop -a 'LINESTRING(10 10, 20 20)' -b 'POLYGON((110 10, 120 10, 120
> > > 20, 110 20, 110 10))' intersection
> > > LINESTRING EMPTY
> > > ```
> > >
> > > But hopefully, it shouldn't be too hard to add a special mode to
> > > RectangleIntersection to return the expected geometry.
> > >
> > > + The resulting geometries gets sometimes rearranged differently. See
> for
> > > example unit test 4 of capi::GEOSIntersection. 2 linestrings are
> expected.
> > > They would not appear with rectangle intersection. By the way, I am
> pretty
> > > novice with geometry computation, and I don't really understand why the
> > > linestrings are expected in that case.
> > >
> > > Another example: where a linestring self-intersects, and the result
> will
> > > be split at the intersection point. For some reason, it does not seem
> to
> > > happen when the linestring is fully within the rectangle boundary (I'm
> also
> > > confused why that is the case).
> > >
> > > ```
> > > % geosop -a 'LINESTRING(13 12, 50 50, 50 15, 10 15)' -b 'POLYGON((10
> 10,
> > > 20 10, 20 20, 10 20, 10 10))' intersection
> > > MULTILINESTRING ((13 12, 15.921052631578947 15), (15.921052631578947
> 15,
> > > 20 19.18918918918919), (20 15, 15.921052631578947 15),
> (15.921052631578947
> > > 15, 10 15))
> > > % geosop -a 'LINESTRING(13 12, 50 50, 50 15, 10 15)' -b 'POLYGON((10
> 10,
> > > 20 10, 20 20, 10 20, 10 10))' clipRect
> > > MULTILINESTRING ((13 12, 20 19.189189189189186), (20 15, 10 15))
> > > ```
> > >
> > > Are there other differences between the current intersection operation,
> > > and the rectangle intersection operation? Or any other issue or
> reasons why
> > > that optimization is not enabled?
> > >
> > > Thank you
> > >
> > > Arno
> > >
> > > PS: I attach the patch I used to enable Rectangle Intersection
> > >
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/geos-devel/attachments/20241127/b85abb4d/attachment.htm>