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 <a...@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
> >
