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
diff --git SRC src/geom/Geometry.cpp DST src/geom/Geometry.cpp
index 433933587..8f960c0fd 100644
--- SRC src/geom/Geometry.cpp
+++ DST src/geom/Geometry.cpp
@@ -73,7 +73,7 @@
#endif
#define SHORTCIRCUIT_PREDICATES 1
-//#define USE_RECTANGLE_INTERSECTION 1
+#define USE_RECTANGLE_INTERSECTION 1
#define USE_RELATENG 1
using namespace geos::algorithm;
@@ -603,13 +603,13 @@ Geometry::intersection(const Geometry* other) const
const Envelope* env = getEnvelopeInternal();
Rectangle rect(env->getMinX(), env->getMinY(),
env->getMaxX(), env->getMaxY());
- return RectangleIntersection::clip(*other, rect).release();
+ return RectangleIntersection::clip(*other, rect);
}
if(other->isRectangle()) {
const Envelope* env = other->getEnvelopeInternal();
Rectangle rect(env->getMinX(), env->getMinY(),
env->getMaxX(), env->getMaxY());
- return RectangleIntersection::clip(*this, rect).release();
+ return RectangleIntersection::clip(*this, rect);
}
#endif
