英文:
Geotools splitting raster and merging it again doesn't work
问题
以下是代码部分的中文翻译：
public ArrayList<GridCoverage2D> createTiles(GridCoverage2D grid) throws FactoryException, TransformException {// 代码内容...}
请注意，这只是代码的翻译，不包括代码之外的任何内容。如果您需要更多帮助或有其他问题，请随时提出。
英文:
first time posting a question here so I'm sorry for any error.
We are trying to let the user choose a raster image (.tiff) so that it can be uploaded to hbase for processing and other operations. Since single images can be quite large, we are trying to split them so that the upload process can be hasten, and then we merge them together again to get the original single image.<br><br>From what I understood, the error seems to lie in the split method.
At the moment we are just trying to split the image and rebuild it locally.
This is the code for the split method
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;}
It kinda works because we get a number of tiffs, but we also get some of them that are not being visualized with tools like qgis. The merge method gives us back a tiff image, but then again it is not visualized in qgis.
Tell me if you need more informations.
Thanks!
EDIT
This is the gdalinfo for the original and correct tif
Driver: GTiff/GeoTIFFFiles: C:\Users\Consulente\Desktop\appoggio_rilascio\originale.tifC:\Users\Consulente\Desktop\appoggio_rilascio\originale.tif.aux.xmlSize is 266, 269Coordinate 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,2Origin = (293436.000000000000000,4647354.000000000000000)Pixel Size = (3.000000000000000,-3.000000000000000)Metadata:AREA_OR_POINT=AreaImage Structure Metadata:COMPRESSION=LZWINTERLEAVE=PIXELCorner 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=RedMin=4230.000 Max=18752.000Minimum=4230.000, Maximum=18752.000, Mean=6493.946, StdDev=1353.641NoData Value=0Metadata:STATISTICS_MAXIMUM=18752STATISTICS_MEAN=6493.9456647057STATISTICS_MINIMUM=4230STATISTICS_STDDEV=1353.6413052408STATISTICS_VALID_PERCENT=94.55Band 2 Block=266x3 Type=UInt16, ColorInterp=GreenMin=3979.000 Max=17611.000Minimum=3979.000, Maximum=17611.000, Mean=5971.019, StdDev=1344.279NoData Value=0Metadata:STATISTICS_MAXIMUM=17611STATISTICS_MEAN=5971.0194223549STATISTICS_MINIMUM=3979STATISTICS_STDDEV=1344.2792863221STATISTICS_VALID_PERCENT=94.55Band 3 Block=266x3 Type=UInt16, ColorInterp=BlueMin=2619.000 Max=16659.000Minimum=2619.000, Maximum=16659.000, Mean=4982.758, StdDev=1525.882NoData Value=0Metadata:STATISTICS_MAXIMUM=16659STATISTICS_MEAN=4982.7577379017STATISTICS_MINIMUM=2619STATISTICS_STDDEV=1525.8817457971STATISTICS_VALID_PERCENT=94.55Band 4 Block=266x3 Type=UInt16, ColorInterp=UndefinedNoData Value=0
and this is the gdalinfo from the uncorrect merged tiff (which should have been equal to the original one)
Driver: GTiff/GeoTIFFFiles: C:\Users\Consulente\Desktop\appoggio_rilascio\merged_tif\merged.tiffSize is 266, 269Coordinate 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,1GeoTransform =41.95133415197179, 0, -2.77698308382438e-0512.50778763370388, 3.722213354286106e-05, 0Metadata:AREA_OR_POINT=AreaTIFFTAG_RESOLUTIONUNIT=1 (unitless)TIFFTAG_XRESOLUTION=1TIFFTAG_YRESOLUTION=1Image Structure Metadata:INTERLEAVE=PIXELCorner 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=RedNoData Value=0Band 2 Block=266x8 Type=UInt16, ColorInterp=GreenNoData Value=0Band 3 Block=266x8 Type=UInt16, ColorInterp=BlueNoData Value=0Band 4 Block=266x8 Type=UInt16, ColorInterp=AlphaNoData Value=0
答案1
得分: 1
从查看GdalInfo输出，您已将CRS更改为EPSG:4326并颠倒了纬度/经度值 - 因此，在您发布的代码中可能存在某种方式来写出瓦片的问题。
如果我要尝试这样做，我会查看图像切片教程，其中关键的内部循环如下：
// 遍历我们的瓦片计数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);// 创建瓦片的包围框Envelope envelope = getTileEnvelope(coverageMinX, coverageMinY, geographicTileWidth, geographicTileHeight,targetCRS, i, j);GridCoverage2D finalCoverage = cropCoverage(gridCoverage, envelope);if (this.getTileScale() != null) {finalCoverage = scaleCoverage(finalCoverage);}// 使用AbstractGridFormat的写入器来写出瓦片fileExtension = "png";File tileFile = new File(tileDirectory, i + "_" + j + "." + fileExtension);final WorldImageWriter wiWriter = new WorldImageWriter(tileFile);// 写入PNG的参数final Format oFormat = wiWriter.getFormat();((AbstractGridFormat) oFormat).getWriter(tileFile).write(finalCoverage, null);}}
以及使用CoverageProcessor类实现的cropCoverage方法：
private GridCoverage2D cropCoverage(GridCoverage2D gridCoverage, Envelope envelope) {CoverageProcessor processor = CoverageProcessor.getInstance();// 手动创建我们想要的操作和参数的示例final ParameterValueGroup param = processor.getOperation("CoverageCrop").getParameters();param.parameter("Source").setValue(gridCoverage);param.parameter("Envelope").setValue(envelope);return (GridCoverage2D) processor.doOperation(param);}
请注意，这些代码片段是Java代码，用于处理地理信息数据的图像切片。
英文:
From looking at the GdalInfo output you have changed the CRS to EPSG:4326 and flipped the lat/lon values - so something is going on in the way you write out the tiles which you haven't posted the code for.
If I was trying to do this I would look at the Image Tiling tutorial, the key inner loop is this:
// iterate over our tile countsfor (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 tileEnvelope 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 tilefileExtension = "png";File tileFile = new File(tileDirectory, i + "_" + j + "." + fileExtension);final WorldImageWriter wiWriter = new WorldImageWriter(tileFile);// writing parameters for pngfinal Format oFormat = wiWriter.getFormat();((AbstractGridFormat) oFormat).getWriter(tileFile).write(finalCoverage, null);}}
and cropCoverage implemented using the CoverageProcessor class:
private GridCoverage2D cropCoverage(GridCoverage2D gridCoverage, Envelope envelope) {CoverageProcessor processor = CoverageProcessor.getInstance();// An example of manually creating the operation and parameters we wantfinal ParameterValueGroup param = processor.getOperation("CoverageCrop").getParameters();param.parameter("Source").setValue(gridCoverage);param.parameter("Envelope").setValue(envelope);return (GridCoverage2D) processor.doOperation(param);}
通过集体智慧和协作来改善编程学习和解决问题的方式。致力于成为全球开发者共同参与的知识库，让每个人都能够通过互相帮助和分享经验来进步。
评论