Geotools拆分栅格并再次合并它不起作用

Question

第一次在这里发问题，如有错误，请见谅。

我们正在尝试让用户选择一个光栅图像（.tiff），以便可以将其上传到 hbase 进行处理和其他操作。由于单个图像可能很大，我们正在尝试将它们拆分以便加快上传过程，然后我们再次将它们合并在一起以获得原始的单个图像。

据我了解，错误似乎在于拆分方法。 目前我们只是尝试分割图像并在本地重建它。

这是拆分方法的代码

public ArrayList<GridCoverage2D> createTiles(GridCoverage2D grid) throws FactoryException, TransformException {
    Geometry geom = buildGeometryFromGrid(grid);
    org.locationtech.jts.geom.Envelope env = geom.getEnvelopeInternal();
    Coverage coverage = GeoHash.coverBoundingBoxMaxHashes(env.getMaxY(), env.getMinX(), env.getMinY(), env.getMaxX(), 1);
    ArrayList<String> geohashes = new ArrayList<>();
    String hash = null;
    if (coverage == null) hash = "";
    else {
        for (String singleHash : coverage.getHashes()) {
            hash = singleHash;
            break;
        }
    }
    if (hash.length() == Constants.defaultPrecision) {
        coverage = GeoHash.coverBoundingBox(env.getMaxY(), env.getMinX(), env.getMinY(), env.getMaxX(), Constants.defaultPrecision + 1);
        geohashes = new ArrayList<>(coverage.getHashes());
    }
    else if (hash.length() < Constants.defaultPrecision) {
        coverage = GeoHash.coverBoundingBox(env.getMaxY(), env.getMinX(), env.getMinY(), env.getMaxX(), Constants.defaultPrecision);
        geohashes = new ArrayList<>(coverage.getHashes());
    }
    else {
        if (hash.length() > Constants.spaceLength) geohashes.add(hash.substring(0, hash.length() - 1));
        else if (hash.length() == Constants.spaceLength) geohashes.add(hash);
        else {
            coverage = GeoHash.coverBoundingBox(env.getMaxY(), env.getMinX(), env.getMinY(), env.getMaxX(), hash.length() + 1);
            geohashes = new ArrayList<>(coverage.getHashes());
        }
    }
    System.out.println(geohashes.size());
    List<ReferencedEnvelope> envelopes = calcuateCuts(geohashes);
    ArrayList<GridCoverage2D> tiles = new ArrayList<>();
    for (ReferencedEnvelope envelope : envelopes) {
        tiles.addAll(getTiles(envelope, grid));
    }
    System.out.println(tiles.size());
    return tiles;
}

它有点工作，因为我们得到了许多 tiff，但我们也得到了一些没有被像 qgis 工具可视化的它们。合并方法返回一个 tiff 图像，但在 qgis 中又没有显示出来。

如果您需要更多信息，请告诉我。 谢谢！

编辑

这是原始且正确的 tif

的 gdalinfo

    Driver: GTiff/GeoTIFF
Files: C:\Users\Consulente\Desktop\appoggio_rilascio\originale.tif
       C:\Users\Consulente\Desktop\appoggio_rilascio\originale.tif.aux.xml
Size is 266, 269
Coordinate System is:
PROJCRS["WGS 84 / UTM zone 33N",
    BASEGEOGCRS["WGS 84",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4326]],
    CONVERSION["UTM zone 33N",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",15,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.9996,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",500000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["unknown"],
        AREA["World - N hemisphere - 12┬░E to 18┬░E - by country"],
        BBOX[0,12,84,18]],
    ID["EPSG",32633]]
Data axis to CRS axis mapping: 1,2
Origin = (293436.000000000000000,4647354.000000000000000)
Pixel Size = (3.000000000000000,-3.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (  293436.000, 4647354.000) ( 12d30'28.04"E, 41d57' 4.05"N)
Lower Left  (  293436.000, 4646547.000) ( 12d30'29.05"E, 41d56'37.91"N)
Upper Right (  294234.000, 4647354.000) ( 12d31' 2.66"E, 41d57' 4.80"N)
Lower Right (  294234.000, 4646547.000) ( 12d31' 3.68"E, 41d56'38.66"N)
Center      (  293835.000, 4646950.500) ( 12d30'45.86"E, 41d56'51.36"N)
Band 1 Block=266x3 Type=UInt16, ColorInterp=Red
  Min=4230.000 Max=18752.000
  Minimum=4230.000, Maximum=18752.000, Mean=6493.946, StdDev=1353.641
  NoData Value=0
  Metadata:
    STATISTICS_MAXIMUM=18752
    STATISTICS_MEAN=6493.9456647057
    STATISTICS_MINIMUM=4230
    STATISTICS_STDDEV=1353.6413052408
    STATISTICS_VALID_PERCENT=94.55
Band 2 Block=266x3 Type=UInt16, ColorInterp=Green
  Min=3979.000 Max=17611.000
  Minimum=3979.000, Maximum=17611.000, Mean=5971.019, StdDev=1344.279
  NoData Value=0
  Metadata:
    STATISTICS_MAXIMUM=17611
    STATISTICS_MEAN=5971.0194223549
    STATISTICS_MINIMUM=3979
    STATISTICS_STDDEV=1344.2792863221
    STATISTICS_VALID_PERCENT=94.55
Band 3 Block=266x3 Type=UInt16, ColorInterp=Blue
  Min=2619.000 Max=16659.000
  Minimum=2619.000, Maximum=16659.000, Mean=4982.758, StdDev=1525.882
  NoData Value=0
  Metadata:
    STATISTICS_MAXIMUM=16659
    STATISTICS_MEAN=4982.7577379017
    STATISTICS_MINIMUM=2619
    STATISTICS_STDDEV=1525.8817457971
    STATISTICS_VALID_PERCENT=94.55
Band 4 Block=266x3 Type=UInt16, ColorInterp=Undefined
  NoData Value=0

这是来自不正确合并 tiff 的gdalinfo（应该与原始tiff 相同）

    Driver: GTiff/GeoTIFF
Files: C:\Users\Consulente\Desktop\appoggio_rilascio\merged_tif\merged.tiff
Size is 266, 269
Coordinate System is:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4326]]
Data axis to CRS axis mapping: 2,1
GeoTransform =
  41.95133415197179, 0, -2.77698308382438e-05
  12.50778763370388, 3.722213354286106e-05, 0
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_RESOLUTIONUNIT=1 (unitless)
  TIFFTAG_XRESOLUTION=1
  TIFFTAG_YRESOLUTION=1
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (  41.9513342,  12.5077876) ( 41d57' 4.80"E, 12d30'28.04"N)
Lower Left  (  41.9438641,  12.5077876) ( 41d56'37.91"E, 12d30'28.04"N)
Upper Right (  41.9513342,  12.5176887) ( 41d57' 4.80"E, 12d31' 3.68"N)
Lower Right (  41.9438641,  12.5176887) ( 41d56'37.91"E, 12d31' 3.68"N)
Center      (  41.9475991,  12.5127382) ( 41d56'51.36"E, 12d30'45.86"N)
Band 1 Block=266x8 Type=UInt16, ColorInterp=Red
  NoData Value=0
Band 2 Block=266x8 Type=UInt16, ColorInterp=Green
  NoData Value=0
Band 3 Block=266x8 Type=UInt16, ColorInterp=Blue
  NoData Value=0
Band 4 Block=266x8 Type=UInt16, ColorInterp=Alpha
  NoData Value=0

Answer 1

通过查看 GdalInfo 输出，您已将 CRS 更改为 EPSG:4326 并翻转了纬度/经度值 - 所以在您写出尚未发布代码的图块的方式上发生了一些事情。

如果我想这样做，我会查看Image Tiling tutorial，关键的内部循环是这样的：

// iterate over our tile counts
for (int i = 0; i < htc; i++) {
  for (int j = 0; j < vtc; j++) {

    System.out.println("Processing tile at indices i: " + i + " and j: " + j);
    // create the envelope of the tile
    Envelope envelope = getTileEnvelope(coverageMinX, coverageMinY, geographicTileWidth, geographicTileHeight,
        targetCRS, i, j);

    GridCoverage2D finalCoverage = cropCoverage(gridCoverage, envelope);

    if (this.getTileScale() != null) {
      finalCoverage = scaleCoverage(finalCoverage);
    }

    // use the AbstractGridFormat's writer to write out the tile
    fileExtension = "png";
    File tileFile = new File(tileDirectory, i + "_" + j + "." + fileExtension);
    final WorldImageWriter wiWriter = new WorldImageWriter(tileFile);

    // writing parameters for png
    final Format oFormat = wiWriter.getFormat();
    ((AbstractGridFormat) oFormat).getWriter(tileFile).write(finalCoverage, null);
  }
}

和cropCoverage 使用CoverageProcessor 类实现：

private GridCoverage2D cropCoverage(GridCoverage2D gridCoverage, Envelope envelope) {
    CoverageProcessor processor = CoverageProcessor.getInstance();

    // An example of manually creating the operation and parameters we want
    final ParameterValueGroup param = processor.getOperation("CoverageCrop").getParameters();
    param.parameter("Source").setValue(gridCoverage);
    param.parameter("Envelope").setValue(envelope);

    return (GridCoverage2D) processor.doOperation(param);
  }