Output from gdal_translate projwin is different from geoserver WCS GetCoverage request?

Question

I've a hard time understanding which one should be the correct result by gdal_translate or GeoServer.

My system (Ubuntu 22.04) has:

• GDAL 3.4.1
• Geoserver 2.22 running on Tomcat on port 9090

I have a small tiff file (link to download at the bottom of this thread) in EPSG:4326, then I wanted to have a small subset on Lat and Long axes from it by running:

 gdal_translate -projwin 137.915315044103 -36.51629558851893 154.64770004642423 -40.081875882617666 test.tif gdal.tif

 curl 'http://localhost:9090/geoserver/wcs?service=WCS&version=2.0.1&request=GetCoverage&format=image/tiff&coverageId=test_raster:output&subset=Lat(-40.081875882617666,-36.51629558851893)&subset=Long(137.915315044103,154.64770004642423)' -o geoserver.tiff

• Then gdalinfo gdal.tif returns:

Size is 17, 4

and gdalinfo geoserver.tif returns

Size is 17, 3

• Also the geo bbox of gdal.tif and geoserver.tif are different. First file gives:

Upper Left  ( 136.9750000, -35.9750000)

Second file gives:

Upper Left  ( 137.9750000, -36.9750000)

More importantly, the pixels are shifted from gdal.tif at the top to geoserver.tif at the bottom which you can see from this image:

Here is the link to download test.tif, gdal.tif and geoserver.tif https://drive.google.com/file/d/1dawCy38FMwx4J7OQsky-gAHHOXsMuDcy/view?usp=sharing

Which result (gdal.tif or geoserver.tif) should be correct?

Full gdalinfo for gdal.tif

Driver: GTiff/GeoTIFF
Files: gdal.tif
Size is 17, 4
Coordinate System is:
GEOGCRS["WGS 84",
   Data axis to CRS axis mapping: 2,1
Origin = (136.974999999999994,-35.975000000000001)
Pixel Size = (1.000000000000000,-1.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  ( 136.9750000, -35.9750000) (136d58'30.00"E, 35d58'30.00"S)
Lower Left  ( 136.9750000, -39.9750000) (136d58'30.00"E, 39d58'30.00"S)
Upper Right ( 153.9750000, -35.9750000) (153d58'30.00"E, 35d58'30.00"S)
Lower Right ( 153.9750000, -39.9750000) (153d58'30.00"E, 39d58'30.00"S)
Center      ( 145.4750000, -37.9750000) (145d28'30.00"E, 37d58'30.00"S)
Band 1 Block=17x4 Type=Byte, ColorInterp=Gray
  NoData Value=0

and for geoserver.tif:

Driver: GTiff/GeoTIFF
Files: geoserver.tif
Size is 17, 3
Coordinate System is:
GEOGCRS["WGS 84",
   Data axis to CRS axis mapping: 2,1
Origin = (137.974999999999994,-36.975000000000001)
Pixel Size = (1.000000000000000,-1.000000000000000)
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_RESOLUTIONUNIT=1 (unitless)
  TIFFTAG_XRESOLUTION=1
  TIFFTAG_YRESOLUTION=1
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  ( 137.9750000, -36.9750000) (137d58'30.00"E, 36d58'30.00"S)
Lower Left  ( 137.9750000, -39.9750000) (137d58'30.00"E, 39d58'30.00"S)
Upper Right ( 154.9750000, -36.9750000) (154d58'30.00"E, 36d58'30.00"S)
Lower Right ( 154.9750000, -39.9750000) (154d58'30.00"E, 39d58'30.00"S)
Center      ( 146.4750000, -38.4750000) (146d28'30.00"E, 38d28'30.00"S)
Band 1 Block=16x16 Type=Byte, ColorInterp=Gray
  NoData Value=0

UPDATE 2: I tested with geotools 29 snapshot with this code below and it returns grid size 18 x 5 instead

        File file = new File("/home/aaa/Downloads/temp/test.tif");

        Hints hint = new Hints();
        hint.put(Hints.DEFAULT_COORDINATE_REFERENCE_SYSTEM, WGS84);
        hint.put(Hints.FORCE_LONGITUDE_FIRST_AXIS_ORDER, Boolean.TRUE);

        GeoTiffReader reader = new GeoTiffReader(file, hint);
        GridCoverage2D coverage = reader.read(null);

        CoverageProcessor processor = CoverageProcessor.getInstance();

        final ParameterValueGroup param = processor.getOperation("CoverageCrop").getParameters();

        // 137.915315044103 -36.51629558851893 154.64770004642423 -40.081875882617666
        double minLong = 137.915315044103;
        double maxLong = 154.64770004642423;
        double minLat = -40.081875882617666;
        double maxLat = -36.51629558851893;

        final GeneralEnvelope crop = new GeneralEnvelope(new GeographicBoundingBoxImpl(minLong, maxLong,
                minLat, maxLat));

        ReferencedEnvelope envelope = new ReferencedEnvelope(new Envelope(minLong, maxLong,
                minLat, maxLat), WGS84);

        param.parameter("Source").setValue(coverage);
        param.parameter("Envelope").setValue(envelope);

        GridCoverage2D cropped = (GridCoverage2D) processor.doOperation(param);

        GridCoverageFactory gcf = new GridCoverageFactory();
        GridCoverage2D gc = gcf.create("name", cropped.getRenderedImage(), cropped.getEnvelope());
        String url = "/tmp/geotools.tif";
        File outputFile = new File(url);
        GeoTiffWriter writer = new GeoTiffWriter(outputFile);
        writer.write(gc, null);
        writer.dispose();

with gdalinfo

gdalinfo /tmp/geotools.tif 
Size is 18, 5
Coordinate System is:
GEOGCRS["WGS 84",
Data axis to CRS axis mapping: 2,1
Origin = (136.974999999999994,-35.975000000000001)
Pixel Size = (1.000000000000000,-1.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Corner Coordinates:
Upper Left  ( 136.9750000, -35.9750000) (136d58'30.00"E, 35d58'30.00"S)
Lower Left  ( 136.9750000, -40.9750000) (136d58'30.00"E, 40d58'30.00"S)
Upper Right ( 154.9750000, -35.9750000) (154d58'30.00"E, 35d58'30.00"S)
Lower Right ( 154.9750000, -40.9750000) (154d58'30.00"E, 40d58'30.00"S)
Center      ( 145.9750000, -38.4750000) (145d58'30.00"E, 38d28'30.00"S)
Band 1 Block=18x5 Type=Byte, ColorInterp=Gray
  NoData Value=0

Answer 1

There is no shift in the pixels of the GDAL TIFF and the GeoServer TIFF if the images are looked as georeferenced for example with QGIS. The difference between 3 and 4 rows comes from the way how the programs deal with the pixels which do not fall completely inside the projwin/subset. GDAL includes the intersected pixels into the result while GeoServer WCS does not.

The behavior of the -projwin in gdal_translate is defined by the GDAL developers so probably it is "correct", even there has been some changes during the lifetile

In GDAL 2.1.0 and 2.1.1, using -projwin with coordinates not aligned with pixels will result in a sub-pixel shift. This has been corrected in later versions. When selecting non-nearest neighbour resampling, starting with GDAL 2.1.0, sub-pixel accuracy is however used to get better results.

GeoServer WCS should act as is defined in the WCS standard.

A.1.32 GetCoverage trimming within coverage limits Test id: /conf/core/getCoverage-request-trim-within-extent

Test Purpose: Requirement /req/core/getCoverage-request-trim-within-extent:

Let the extent of the coverage’s gml:Envelope along the dimension specified in the trim request range from L to H. Then, for the trim bounds trimLow and trimHigh the following shall hold:

L <= trimLow <= trimHigh <= H.

In the GeoServer response the trimHigh value is larger than the maximum latitude

trimHigh = -36.51629558851893
H        = -36.9750000 

It looks that the GeoServer response is not correct by the WCS 2.0 standard. Comparison with gdal_translate that is using the GTiff driver is not really relevant.

Answer 2

I considered the result from geoserver is the correct one and I confirmed that with R raster package.

Fortunately, I found the formula from R raster which I could understand. Here is the R code to return the same subset like geoserver WCS GetCoverage subset.

Copied from:

• https://github.com/rspatial/raster/blob/f71e83b9838ca79c05c4b0a881b970966b2d1191/R/xyCell.R#L95
• https://github.com/rspatial/raster/blob/805e5e80f0a7363e583b2486c914bdd300b47728/R/crop.R#L53

suppressPackageStartupMessages(library(raster))

x <- raster(ncols=44, nrows=36)

# min_long, min_lat, max_long, max_lat

xmin <- 111.975
ymin <- -44.975
xmax <- 155.975
ymax <- -8.975

extent(x) <- extent(xmin, xmax, ymin, ymax)

# min_long, min_lat, max_long, max_lat

input_xmin <- {}
input_ymin <- {}
input_xmax <- {}
input_ymax <- {}

out <- extent(input_xmin, input_xmax, input_ymin, input_ymax)

col1 <- colFromX(x, xmin(out)+0.5*xres(x))
col2 <- colFromX(x, xmax(out)-0.5*xres(x))


# gridIndex1 < gridIndex2

row1 <- rowFromY(x, ymax(out)-0.5*yres(x))
row2 <- rowFromY(x, ymin(out)+0.5*yres(x))

col1 <- col1 - 1
col2 <- col2 - 1
row1 <- row1 - 1
row2 <- row2 - 1

pixel_x <- col2 - col1 + 1
pixel_y <- row2 - row1 + 1

geo_lower_bound_x <- xmin + (col1) * xres(x)
geo_lower_bound_y <- ymax - (row2 + 1) * yres(x)
geo_upper_bound_x <- xmin + (col2 + 1) * xres(x)
geo_upper_bound_y <- ymax - (row1) * yres(x)

cat(pixel_x, pixel_y, geo_lower_bound_x, geo_lower_bound_y, geo_upper_bound_x, geo_upper_bound_y)

And the python test code to calculate the grid extent from the geo-referenced extent based on the above formula:

import os

# nrow = 10, col = 10, ncell = 100
# resx = 36, resy = 18
# extent = -180, 180, -90, 90 (xmin, xmax, ymin, ymax)

xmin = -180
xmax = 180
ymin = -90
ymax = 90
xres = 36
yres = 18

def colFromX(x):
    colnr = int( (x - xmin) / xres ) + 1
    return colnr

def rowFromY(y):
    rownr = 1 + ( int( (ymax - y) / yres ) )
    return rownr
    

def cal(subset_xmin, subset_ymin, subset_xmax, subset_ymax):
    col1 = colFromX(subset_xmin + 0.5 * xres) 
    col2 = colFromX(subset_xmax - 0.5 * xres)
    row1 = rowFromY(subset_ymax - 0.5 * yres) 
    row2 = rowFromY(subset_ymin + 0.5 * yres)
    
    print(col1 - 1, col2 - 1, row1 - 1, row2 - 1)