H3空间索引的使用（Java）

首先添加Maven依赖：

<dependencies>
    <!-- H3 Core Library -->
    <dependency>
        <groupId>com.uber.h3</groupId>
        <artifactId>h3</artifactId>
        <version>3.7.2</version>
    </dependency>
    
    <!-- GeoTools -->
    <dependency>
        <groupId>org.geotools</groupId>
        <artifactId>gt-shapefile</artifactId>
        <version>25.2</version>
    </dependency>
    <dependency>
        <groupId>org.geotools</groupId>
        <artifactId>gt-main</artifactId>
        <version>25.2</version>
    </dependency>
</dependencies>

Java实现代码：

import com.uber.h3core.H3Core;
import org.geotools.data.FileDataStore;
import org.geotools.data.FileDataStoreFinder;
import org.geotools.data.simple.SimpleFeatureCollection;
import org.geotools.data.simple.SimpleFeatureIterator;
import org.geotools.data.simple.SimpleFeatureSource;
import org.locationtech.jts.geom.Coordinate;
import org.locationtech.jts.geom.Geometry;
import org.locationtech.jts.geom.Point;
import org.opengis.feature.simple.SimpleFeature;

import java.io.File;
import java.io.IOException;
import java.util.*;

public class H3SpatialIndex {

    private final H3Core h3;
    private final Map<String, List<PointData>> index = new HashMap<>();
    private final int resolution = 9; // 调整分辨率（0-15）

    public H3SpatialIndex() throws IOException {
        h3 = H3Core.newInstance();
    }

    // 加载Shapefile数据并建立索引
    public void buildIndex(String shpPath) throws Exception {
        File file = new File(shpPath);
        FileDataStore store = FileDataStoreFinder.getDataStore(file);
        SimpleFeatureSource featureSource = store.getFeatureSource();
        SimpleFeatureCollection features = featureSource.getFeatures();

        try (SimpleFeatureIterator iterator = features.features()) {
            while (iterator.hasNext()) {
                SimpleFeature feature = iterator.next();
                Geometry geom = (Geometry) feature.getDefaultGeometry();
                
                if (geom instanceof Point) {
                    Point point = (Point) geom;
                    double lat = point.getY();
                    double lng = point.getX();
                    
                    // 生成H3索引
                    String h3Address = h3.geoToH3Address(lat, lng, resolution);
                    
                    // 存储数据
                    PointData data = new PointData(
                        feature.getID(),
                        lat,
                        lng,
                        feature.getAttributes()
                    );
                    
                    index.computeIfAbsent(h3Address, k -> new ArrayList<>()).add(data);
                }
            }
        }
        store.dispose();
    }

    // 半径查询（单位：米）
    public List<PointData> query(double centerLat, double centerLng, double radius) {
        List<PointData> results = new ArrayList<>();
        
        // 确定搜索半径对应的k值
        int k = estimateK(radius);
        
        // 获取中心点H3地址
        String centerH3 = h3.geoToH3Address(centerLat, centerLng, resolution);
        
        // 获取k-ring所有相邻单元
        List<String> neighbors = h3.kRing(centerH3, k);
        
        // 收集所有候选点
        for (String h3Address : neighbors) {
            List<PointData> points = index.getOrDefault(h3Address, Collections.emptyList());
            for (PointData p : points) {
                if (haversineDistance(centerLat, centerLng, p.lat, p.lng) <= radius) {
                    results.add(p);
                }
            }
        }
        
        return results;
    }

    // 估算k值（简单实现，实际需根据H3分辨率调整）
    private int estimateK(double radius) {
        // 分辨率对应的六边形边长参考：
        // https://h3geo.org/docs/core-library/restable
        // 例如分辨率9对应边长约174米
        return (int) Math.ceil(radius / 174);
    }

    // Haversine距离计算
    private static double haversineDistance(double lat1, double lng1, double lat2, double lng2) {
        final int R = 6371; // 地球半径（千米）
        double dLat = Math.toRadians(lat2 - lat1);
        double dLng = Math.toRadians(lng2 - lng1);
        double a = Math.sin(dLat/2) * Math.sin(dLat/2) +
                   Math.cos(Math.toRadians(lat1)) * Math.cos(Math.toRadians(lat2)) *
                   Math.sin(dLng/2) * Math.sin(dLng/2);
        double c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1-a));
        return R * c * 1000; // 转换为米
    }

    public static class PointData {
        public final String id;
        public final double lat;
        public final double lng;
        public final List<Object> attributes;

        public PointData(String id, double lat, double lng, List<Object> attributes) {
            this.id = id;
            this.lat = lat;
            this.lng = lng;
            this.attributes = attributes;
        }
    }

    public static void main(String[] args) throws Exception {
        H3SpatialIndex spatialIndex = new H3SpatialIndex();
        
        // 1. 构建索引
        spatialIndex.buildIndex("path/to/your/points.shp");
        
        // 2. 示例查询（天安门坐标）
        double centerLat = 39.9087;
        double centerLng = 116.3975;
        double radius = 1000; // 1公里
        
        List<PointData> results = spatialIndex.query(centerLat, centerLng, radius);
        
        // 输出结果
        System.out.println("找到 " + results.size() + " 个点：");
        for (PointData p : results) {
            System.out.printf("ID: %s, 坐标: (%.6f, %.6f)\n", p.id, p.lat, p.lng);
        }
    }
}

使用说明：

1. 需要根据实际情况调整以下参数：

   • resolution：H3分辨率（0-15），值越大网格越小

   • estimateK()：根据实际分辨率调整k值估算逻辑

   • 确保Shapefile使用WGS84坐标系（EPSG:4326）

2. 功能特点：

   • 使用GeoTools读取Shapefile

   • 建立H3空间索引

   • 支持半径范围查询

   • 使用Haversine公式进行精确距离过滤

3. 性能优化建议：

   • 预处理时存储不同分辨率的索引

   • 使用空间数据库进行持久化存储

   • 对大量数据使用并发处理

如果需要处理非WGS84坐标系的Shapefile，需要添加坐标转换代码：

import org.geotools.referencing.CRS;
import org.opengis.referencing.crs.CoordinateReferenceSystem;
import org.opengis.referencing.operation.MathTransform;

// 在buildIndex方法中添加坐标转换：
CoordinateReferenceSystem sourceCRS = featureSource.getSchema().getCoordinateReferenceSystem();
CoordinateReferenceSystem targetCRS = CRS.decode("EPSG:4326");
MathTransform transform = CRS.findMathTransform(sourceCRS, targetCRS);

Geometry targetGeometry = JTS.transform(geom, transform);
Point point = (Point) targetGeometry;

GeoTools需要配置Maven仓库：

<repositories>
    <repository>
        <id>osgeo</id>
        <name>OSGeo Release Repository</name>
        <url>https://repo.osgeo.org/repository/release/</url>
    </repository>
</repositories>