gdalwarp sum resampling algorithm double counting at some specific resolutions

Question

GDALWarp utility is double-counting my pixels in the sum algorithm for resampling in some cases. I don't believe it's something I'm obviously doing wrong.

I have a file, at 30 arcsecond resolution for the globe.

raster = rasterio.open(f'data.tif')
raster = raster.read(1)
print(np.sum(raster))
>>> 1,028,220,143

Now, I want to resample this to different resolutions.

gdalwarp -co "COMPRESS=DEFLATE" -r sum -dstnodata 0 -tr 0.6333333 0.6333333 -te -180 -90 180 90 data.tif data_0.63deg.tif
>>> Creating output file that is 568P x 284L.

The result I get is:

raster = rasterio.open(f'data_0.63deg.tif')
raster = raster.read(1)
print(np.sum(raster))
>>> 1,028,220,142.9999847

Great, this is what I expect.

gdalwarp -co "COMPRESS=DEFLATE" -r sum -dstnodata 0 -tr 0.6333 0.51 -te -180 -90 180 90 data.tif data_0.63-0.51deg.tif
>>> Creating output file that is 568P x 353L.

The result I get is:

raster = rasterio.open(f'data_0.63-0.51deg.tif')
raster = raster.read(1)
print(np.sum(raster))
>>> 1,028,220,143.0000087

Still great, this is what I expect.

And when I get to some resolutions, the sum resampling algorithm double counts everything somehow.

gdalwarp -co "COMPRESS=DEFLATE" -r sum -dstnodata 0 -tr 0.5 0.5 -te -180 -90 180 90 data.tif data_0.5deg.tif
>>> Creating output file that is 720P x 360L.

The result I get is:

raster = rasterio.open(f'data_0.5deg.tif')
raster = raster.read(1)
print(np.sum(raster))
>>> 2,056,440,285.9999878

So, every pixel is now double counted. I thought it may be a weird cell alignment issue, but trying to do slightly different gives the same doubling:

gdalwarp -co "COMPRESS=DEFLATE" -r sum -dstnodata 0 -tr 0.51 0.51 -te -180 -90 180 90 data.tif data_0.51deg.tif
>>> Creating output file that is 706P x 353L.

The result I get is:

raster = rasterio.open(f'data_0.51deg.tif')
raster = raster.read(1)
print(np.sum(raster))
>>> 2,056,440,285.9999878

My GDAL version is:

$gdalwarp --version
GDAL 3.1.2, released 2020/07/07

My OS is MacOS 10.15.4 Catalina

GDALinfo on my data.tif:

$gdalinfo data.tif
Driver: GTiff/GeoTIFF
Files: data.tif
Size is 43200, 21600
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
Origin = (-180.000000000000000,90.000000000000000)
Pixel Size = (0.008333333000000,-0.008333333000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=DEFLATE
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (-180.0000000,  90.0000000) (180d 0' 0.00"W, 90d 0' 0.00"N)
Lower Left  (-180.0000000, -89.9999928) (180d 0' 0.00"W, 89d59'59.97"S)
Upper Right ( 179.9999856,  90.0000000) (179d59'59.95"E, 90d 0' 0.00"N)
Lower Right ( 179.9999856, -89.9999928) (179d59'59.95"E, 89d59'59.97"S)
Center      (  -0.0000072,   0.0000036) (  0d 0' 0.03"W,  0d 0' 0.01"N)
Band 1 Block=43200x1 Type=Float64, ColorInterp=Gray
  NoData Value=0

I must say that I am at a loss now. Anyone has ever encountered something like that? Any suggestion for getting the same resampling using maybe rasterio alone and bypassing what may be a bug in the gdalwarp utility? Or maybe I'm doing something very obviously wrong that I am completely missing.

Note also that I have tried to do some manual debugging by opening the data.tif (30 arcsecond) file in QGIS, and opening the data_0.5deg.tif too. I identified a location where I had a low enough numbers of 30 arcsecond pixel to read in a 0.5deg cell, and manually summed them. I obtain the same value as what I expected, no doubling. So each pixel individually (I was only able to look at two such 0.5deg cells) is correct, yet the sum double counts somehow. I also tried the sum in QGIS using the "raster layer statistics" process in the toolbox, with the same double counting appearing.

I also thought "maybe something went wrong, let's rerun", and reran everything in a different order, at different resolutions, and still the problem appear for the same resolutions (I have seen it happen in my case for 0.5deg, 0.51deg, and 0.25deg, for example), while it still does not happen for the other resolutions I have tried.

Answer 1

I was finally able to see what is happening, although I have yet to figure out how to solve it.

When lowering the resolution to 0.5 degree, the sum algorithm works great and as expected everywhere, but then, in the very last pixels column (i.e., for the world, the pixels at 179.5 to 180 and -90 to 90), the sum of the latitude (pixels row) is added.

This is what creates my issue, where locally, everything looks as expected, but globally, there is a doubling. I am very uncertain as to what generate this behavior, but will edit this post as soon as I figure it out.