Geotools splitting raster and merging it again doesn't work

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/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

and this is the gdalinfo from the uncorrect merged tiff (which should have been equal to the original one)

    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

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 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);
  }
}

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 want
    final ParameterValueGroup param = processor.getOperation("CoverageCrop").getParameters();
    param.parameter("Source").setValue(gridCoverage);
    param.parameter("Envelope").setValue(envelope);

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