如何用GeoTools实现网格与多边形空间连接？及点在多边形快速匹配

Hey there! I totally get the pain of slow naive point-in-polygon checks—700k points against 260 polygons is no joke, and iterating every point against every polygon is definitely going to drag. Let's look at the optimized approaches GeoTools offers that'll get you performance closer to R's over function.

1. Use PreparedGeometry for Fast Repeated Contains Checks #

GeoTools provides PreparedGeometry specifically to speed up repeated spatial operations (like checking if multiple points are inside a single polygon). It pre-processes the polygon into an indexed structure, making each subsequent contains() check way faster than the raw Geometry method.

Here's a quick example of how to implement this:

import org.geotools.geometry.jts.JTS;
import org.locationtech.jts.geom.Geometry;
import org.locationtech.jts.geom.Point;
import org.locationtech.jts.prep.PreparedGeometry;
import org.locationtech.jts.prep.PreparedGeometryFactory;

// First, pre-process all polygons into PreparedGeometry objects
List<PreparedGeometry> preparedPolygons = new ArrayList<>();
PreparedGeometryFactory factory = new PreparedGeometryFactory();
for (Geometry polygon : yourPolygonList) {
    preparedPolygons.add(factory.create(polygon));
}

// Then, iterate through your points and check against prepared geometries
for (Point point : yourPointList) {
    for (PreparedGeometry prepGeom : preparedPolygons) {
        if (prepGeom.contains(point)) {
            // Assign point to this polygon
            break;
        }
    }
}

This alone can cut down your per-point check time significantly, since the prepared geometry avoids re-computing spatial properties for each check.

2. Add a Spatial Index to Reduce Candidate Polygons #

Even with PreparedGeometry, checking all 260 polygons per point is unnecessary. Use GeoTools/JTS's STRtree (a spatial index) to quickly narrow down which polygons might contain the point, based on their bounding envelopes. This way, you only check a tiny subset of polygons per point.

Combine this with PreparedGeometry for maximum speed:

import org.locationtech.jts.index.strtree.STRtree;
import org.locationtech.jts.geom.Envelope;

// Build the spatial index with your polygons
STRtree spatialIndex = new STRtree();
for (Geometry polygon : yourPolygonList) {
    spatialIndex.insert(polygon.getEnvelopeInternal(), polygon);
}

// Pre-process the candidate polygons into PreparedGeometry once
PreparedGeometryFactory factory = new PreparedGeometryFactory();

// Now process each point
for (Point point : yourPointList) {
    // Get polygons whose envelopes intersect the point's envelope
    List<Geometry> candidatePolygons = spatialIndex.query(point.getEnvelopeInternal());
    // Check only these candidates with PreparedGeometry
    for (Geometry polygon : candidatePolygons) {
        PreparedGeometry prepGeom = factory.create(polygon); // Or cache these prepped geometries too!
        if (prepGeom.contains(point)) {
            // Assign point to this polygon
            break;
        }
    }
}

Pro tip: Cache the PreparedGeometry for each polygon instead of creating it on the fly—this saves even more time.

3. Batch Processing with Feature Collections & Spatial Joins #

If you're working with Feature objects (instead of raw JTS geometries), you can use GeoTools' built-in spatial join utilities, which handle indexing and optimization under the hood. This is closest to how R's over function works, as it abstracts away the low-level indexing.

Here's a high-level outline:

• Convert your points and polygons into SimpleFeatureCollection instances.
• Create a spatial filter that checks "point is within polygon".
• Use DefaultQuery and FeatureSource to perform the join, which GeoTools optimizes with spatial indexes automatically.

This approach is great if you're already working with GeoTools' feature model, as it reduces boilerplate code.

Why This Works (And How It Compares to R's over) #

R's over function relies on spatial indexing (like R-tree) and optimized geometric operations under the hood—exactly what we're implementing here with GeoTools. By combining spatial indexing to reduce candidates and PreparedGeometry to speed up the actual contains check, you'll see a massive performance boost over the naive approach.

For your dataset (700k points, 260 polygons), the spatial index + prepared geometry combo should bring your runtime down from hours to minutes (or even seconds, depending on your hardware).

内容的提问来源于stack exchange，提问作者user2976698