0

我正在尝试使用 Geotools 用 H3 网格定义的瓷砖平铺历史政治领土。我的输入 shapefile 使用 WGS84/EPSG:4326 CRS。H3 六角网格是面向中心的日晷投影，但我假设它们仍然使用 WGS84 CRS。（只是因为这就是 Kevin Sahr 在启发 H3 项目的 DGGS 软件中使用的，据我所知）。

我得到了一些奇怪的结果，我怀疑问题是因为我没有在我的方法中处理指针投影；既不放置六边形也不将生成的多边形保存在 shapefile 输出中。

这是我的方法：

  1. 摄取区域的 shapefile（多边形和多多边形）
  2. h3.polyfillAddress()提取边界使用该方法以分辨率 3（或 4，在某些情况下）收集适合该边界的六边形
  3. 使用 Geotools 为该组六边形构建并保存 shapefileSimpleFeatureTypeBuilder

1815 年 H3 瓦片在 EPSG:4326 (WGS84) 中亚汗国地图上的未解决投影

输入多边形的 geojson 输入位于：https ://pastebin.com/wH1SgCY8

H3 瓦片的 geojson 输出位于：https ://pastebin.com/0Zuc7awv

地图中心附近的形状文件 - 使用 QGIS 查看的本初子午线/赤道 - 似乎是正确的，但离中心越远我得到奇怪的结果。投影显然是一个问题，但我不知道奇怪的 polyfill 结果是由于 H3 的 polyfill 方法的某些限制还是因为在尝试 polyfill 操作之前应该已经处理了输入边界的投影。

我的 SimpleFeatureTypeBuilder 使用 WGS84 CRS 指定输出。在构建每个多边形之前，我是否应该对每个多边形进行从指针投影（到？）的一些转换？

这是代码：

   private static void saveShapefile(String name, String abb, Set<String> hexes) throws IOException {
    /*
    Takes a country name, it's abbreviation and a set of hexagons that will represent it and writes out a shapefile

     */
    String filename = "src/main/data/hexMaps/politicalBoundaryHexes_" + yearString + "_" + name + ".shp";
    File file = new File(filename);

    H3Core h3 = H3Core.newInstance();
    List<SimpleFeature> features = new ArrayList<>();
    GeometryFactory gf = JTSFactoryFinder.getGeometryFactory();

    DefaultFeatureCollection output = new DefaultFeatureCollection();

    final SimpleFeatureType HEX = createFeatureType();
    SimpleFeatureBuilder featureBuilder = new SimpleFeatureBuilder(HEX);

    ShapefileDataStoreFactory dataStoreFactory = new ShapefileDataStoreFactory();
    Map<String, Serializable> params = new HashMap<>();
    params.put("url", file.toURI().toURL());
    params.put("create spatial index", Boolean.TRUE);

    ShapefileDataStore dataStore = (ShapefileDataStore) dataStoreFactory.createNewDataStore(params);
    dataStore.createSchema(HEX);

    int index = 0;
    for (String h:hexes) {
        List<GeoCoord> geocoords = h3.h3ToGeoBoundary(h);
        int numVerts = geocoords.size();
        Coordinate[] coordinates = new Coordinate[numVerts + 1];
        if(debugging){System.out.println("Processing hex " + index + ", @ " + h);}
        for (int i=0; i<geocoords.size(); i++){
            coordinates[i] = new Coordinate(geocoords.get(i).lat, geocoords.get(i).lng);
        }
        coordinates[numVerts] = new Coordinate(geocoords.get(0).lat, geocoords.get(0).lng);
        Polygon polygon = gf.createPolygon(coordinates);
        Object[] values = new Object[]{polygon, index, h};
        index++;
        featureBuilder.addAll(values);

        SimpleFeature polyFeature = featureBuilder.buildFeature(h);
        features.add(polyFeature);
    }

    org.geotools.data.Transaction transaction = new DefaultTransaction("create");
    String typeName = dataStore.getTypeNames()[0];
    SimpleFeatureSource featureSource = dataStore.getFeatureSource(typeName);

    if (featureSource instanceof SimpleFeatureStore) {
        SimpleFeatureStore featureStore = (SimpleFeatureStore) featureSource;
        SimpleFeatureCollection collection = new ListFeatureCollection(HEX, features);
        featureStore.setTransaction(transaction);
        try {
            featureStore.addFeatures(collection);
            transaction.commit();
        } catch (Exception e) {
            e.printStackTrace();
            transaction.rollback();
        } finally {
            transaction.close();
        }
    } else {
        System.out.println("Writing the Shapefile failed.");
        System.exit(1); // Failure!
    }
}

private static SimpleFeatureType createFeatureType() {
    SimpleFeatureTypeBuilder builder = new SimpleFeatureTypeBuilder();
    builder.setName("polygon");
    builder.setCRS(DefaultGeographicCRS.WGS84);
    builder.add("the_geom", Polygon.class);
    builder.add("id", Long.class);
    builder.add("address", String.class);
    final SimpleFeatureType TYPE = builder.buildFeatureType();
    return TYPE;
}
4

0 回答 0