Discrepancies in COG metadata computed from "save_cog_with_dask" compared to GDAL's

[remi-braun]

I have discrepancies in metadata between COGs saved with odc-geo and GDAL (via rioxarray), and I think this is maybe a bug in odc-geo:

• missing LAYOUT=COG in Image Structure Metadata
• discrepancy between the bands colorinterp: Band 1 is ColorInterp=Gray when the two other have ColorInterp=Undefined
• INTERLEAVE=BAND instead of PIXEL (this is maybe it's a feature ?)

If this is all wanted, is there a way to change it to mimic GDAL's behaviour?

---

Computed with odcgeo 4.10.0:
gdalinfo cog_dask.tif

Size is 3700, 1536
Coordinate System is:
PROJCRS["WGS 84 / Pseudo-Mercator",
    BASEGEOGCRS["WGS 84",
        ENSEMBLE["World Geodetic System 1984 ensemble",
            MEMBER["World Geodetic System 1984 (Transit)"],
            MEMBER["World Geodetic System 1984 (G730)"],
            MEMBER["World Geodetic System 1984 (G873)"],
            MEMBER["World Geodetic System 1984 (G1150)"],
            MEMBER["World Geodetic System 1984 (G1674)"],
            MEMBER["World Geodetic System 1984 (G1762)"],
            MEMBER["World Geodetic System 1984 (G2139)"],
            MEMBER["World Geodetic System 1984 (G2296)"],
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]],
            ENSEMBLEACCURACY[2.0]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4326]],
    CONVERSION["Popular Visualisation Pseudo-Mercator",
        METHOD["Popular Visualisation Pseudo Mercator",
            ID["EPSG",1024]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["False easting",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["easting (X)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["northing (Y)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Web mapping and visualisation."],
        AREA["World between 85.06┬░S and 85.06┬░N."],
        BBOX[-85.06,-180,85.06,180]],
    ID["EPSG",3857]]
Data axis to CRS axis mapping: 1,2
Origin = (6166949.000000000000000,-2382654.500000000000000)
Pixel Size = (0.500000000000000,-0.500000000000000)
Metadata:
  TIFFTAG_XRESOLUTION=1
  TIFFTAG_YRESOLUTION=1
  TIFFTAG_RESOLUTIONUNIT=1 (unitless)
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=DEFLATE
  INTERLEAVE=BAND
  PREDICTOR=2
Corner Coordinates:
Upper Left  ( 6166949.000,-2382654.500) ( 55d23'55.12"E, 20d55'21.44"S)
Lower Left  ( 6166949.000,-2383422.500) ( 55d23'55.12"E, 20d55'44.64"S)
Upper Right ( 6168799.000,-2382654.500) ( 55d24'54.95"E, 20d55'21.44"S)
Lower Right ( 6168799.000,-2383422.500) ( 55d24'54.95"E, 20d55'44.64"S)
Center      ( 6167874.000,-2383038.500) ( 55d24'25.04"E, 20d55'33.04"S)
Band 1 Block=2048x1456 Type=Byte, ColorInterp=Gray
  Description = RED
  NoData Value=0
  Overviews: 1850x768, 925x384
Band 2 Block=2048x1456 Type=Byte, ColorInterp=Undefined
  Description = GREEN
  NoData Value=0
  Overviews: 1850x768, 925x384
Band 3 Block=2048x1456 Type=Byte, ColorInterp=Undefined
  Description = BLUE
  NoData Value=0
  Overviews: 1850x768, 925x384

gdalinfo cog_rioxarray.tif

Size is 3700, 1535
Coordinate System is:
PROJCRS["WGS 84 / Pseudo-Mercator",
    BASEGEOGCRS["WGS 84",
        ENSEMBLE["World Geodetic System 1984 ensemble",
            MEMBER["World Geodetic System 1984 (Transit)"],
            MEMBER["World Geodetic System 1984 (G730)"],
            MEMBER["World Geodetic System 1984 (G873)"],
            MEMBER["World Geodetic System 1984 (G1150)"],
            MEMBER["World Geodetic System 1984 (G1674)"],
            MEMBER["World Geodetic System 1984 (G1762)"],
            MEMBER["World Geodetic System 1984 (G2139)"],
            MEMBER["World Geodetic System 1984 (G2296)"],
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]],
            ENSEMBLEACCURACY[2.0]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4326]],
    CONVERSION["Popular Visualisation Pseudo-Mercator",
        METHOD["Popular Visualisation Pseudo Mercator",
            ID["EPSG",1024]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["False easting",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["easting (X)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["northing (Y)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Web mapping and visualisation."],
        AREA["World between 85.06┬░S and 85.06┬░N."],
        BBOX[-85.06,-180,85.06,180]],
    ID["EPSG",3857]]
Data axis to CRS axis mapping: 1,2
Origin = (6166949.000000000000000,-2382654.500000000000000)
Pixel Size = (0.500000000000000,-0.500000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  LAYOUT=COG
  COMPRESSION=DEFLATE
  INTERLEAVE=PIXEL
  PREDICTOR=2
Corner Coordinates:
Upper Left  ( 6166949.000,-2382654.500) ( 55d23'55.12"E, 20d55'21.44"S)
Lower Left  ( 6166949.000,-2383422.000) ( 55d23'55.12"E, 20d55'44.62"S)
Upper Right ( 6168799.000,-2382654.500) ( 55d24'54.95"E, 20d55'21.44"S)
Lower Right ( 6168799.000,-2383422.000) ( 55d24'54.95"E, 20d55'44.62"S)
Center      ( 6167874.000,-2383038.250) ( 55d24'25.04"E, 20d55'33.03"S)
Band 1 Block=512x512 Type=Byte, ColorInterp=Undefined
  Description = RED
  NoData Value=0
  Overviews: 1850x767, 925x383, 462x191
Band 2 Block=512x512 Type=Byte, ColorInterp=Undefined
  Description = GREEN
  NoData Value=0
  Overviews: 1850x767, 925x383, 462x191
Band 3 Block=512x512 Type=Byte, ColorInterp=Undefined
  Description = BLUE
  NoData Value=0
  Overviews: 1850x767, 925x383, 462x191

Here are the files: cogs.zip

[Kirill888]

@remi-braun

missing LAYOUT=COG in Image Structure Metadata

• LAYOUT=COG is some GDAL specific extra marker stored in "ghost area" of the TIFF image header and also implies extra leader/trailer bytes are stored just before and after EVERY compressed tile that seems to be there to enable detection of "file structure corruption by non-COG aware writers?"
  • we don't add those markers, hence no LAYOUT=COG
  • files are COG though, and unlike GDAL produced files should not have any "ghost area"s

INTERLEAVE=BAND instead of PIXEL (this is maybe it's a feature ?)

• yes it is a feature
• Band,Y,X arrays are stores as separate bands (i.e. INTERLEAVE=BAND)
• Y,X,Band arrays are stored as a single multi-channel band (i.e. INTERLEAVE=PIXEL)
• so perform transpose before saving to get RGB image on output

discrepancy between the bands colorinterp: Band 1 is ColorInterp=Gray when the two other have ColorInterp=Undefined

Again this is GDAL specific interpretation of these TIFF tags:

TiffTag 262 PhotometricInterpretation @104 SHORT @116 = MINISBLACK
TiffTag 338 ExtraSamples @324 SHORT[2] @336 = ('UNSPECIFIED', 'UNSPECIFIED')

⬇ tiffinfo cog_dask.tif

=== TIFF directory 0 ===
TIFF Directory at offset 0x10 (16)
  Image Width: 3700 Image Length: 1536
  Tile Width: 2048 Tile Length: 1456
  Resolution: 1, 1 (unitless)
  Bits/Sample: 8
  Compression Scheme: AdobeDeflate
  Photometric Interpretation: min-is-black
  Extra Samples: 2<unspecified, unspecified>
  Samples/Pixel: 3
  Planar Configuration: separate image planes
  Tag 33550: 0.500000,0.500000,0.000000
  Tag 33922: 0.000000,0.000000,0.000000,6166949.000000,-2382654.500000,0.000000
  Tag 34735: 1,1,0,7,1024,0,1,1,1025,0,1,1,1026,34737,25,0,2049,34737,7,25,2054,0,1,9102,3072,0,1,3857,3076,0,1,9001
  Tag 34737: WGS 84 / Pseudo-Mercator|WGS 84|
  GDAL Metadata: <GDALMetadata><Item name="DESCRIPTION" sample="0" role="description">RED</Item><Item name="DESCRIPTION" sample="1" role="description">GREEN</Item><Item name="DESCRIPTION" sample="2" role="description">BLUE</Item></GDALMetadata>
  GDAL NoDataValue: 0
  Predictor: horizontal differencing 2 (0x2)

in TIFF "extra samples" can only be marked as "2 kinds of Alpha channel, or something else"

#define	TIFFTAG_EXTRASAMPLES		338	/* !info about extra samples */
#define	    EXTRASAMPLE_UNSPECIFIED	0	/* !unspecified data */
#define	    EXTRASAMPLE_ASSOCALPHA	1	/* !associated alpha data */
#define	    EXTRASAMPLE_UNASSALPHA	2	/* !unassociated alpha data */

so you can't have Gray for all the bands in TIFF in general, can't generate such file with GDAL either.

[remi-braun]

Thanks a lot for explaining all this! 🙏

As I am a bit out of my depth here, I'll explain the usecase:

I wouldn't have dig in all this if I hadn't an issue with a viewer: the nodata of the COG generated with odc.geo.save_cog_with_dask was displayed in black, although the COG was perfectly fine in QGIS.

When looking for answers to this issue, I was told that LAYOUT=COG was mandatory to make it work correctly, and also that the colorinterp of all bands should be the same.

If LAYOUT=COG isn't mandatory, how could I set ColorInterp=Undefined to all bands?
(I'm assuming nodata is displayed as black because of this: Photometric Interpretation: min-is-black -> set PHOTOMETRIC.RGB ?)

PS: In the end I had to regenerate the COG with GDAL functions to make it work correctly.

[remi-braun]

I digging a bit more and I was wondering in this snippet of _make_empty_cog:

    ax, yaxis = yaxis_from_shape(shape, gbox)
    im_shape = shape_(shape[yaxis : yaxis + 2])
    photometric = PHOTOMETRIC.MINISBLACK
    planarconfig: Optional[PLANARCONFIG] = PLANARCONFIG.SEPARATE
    if ax == "YX":
        nsamples = 1
    elif ax == "YXS":
        nsamples = shape[-1]
        planarconfig = PLANARCONFIG.CONTIG
        if nsamples in (3, 4):
            photometric = PHOTOMETRIC.RGB
    else:
        nsamples = shape[0]
        if nsamples == 1:
            planarconfig = None

why photometric = PHOTOMETRIC.MINISBLACK for ax == "SXY" (this is what happens for me) ?
Shouldn't it be also RGB if nsamples in (3, 4)?

[Kirill888]

@remi-braun before digging any further, please try saving in RGB mode, that's what you need in this case. By RGB mode I mean arrange your data in Y,X,B order before calling save_cog_with_dask. You should then get proper image with photometric = PHOTOMETRIC.RGB.

and yes we probably need to expose more options like photometric= to the user

[Kirill888]

Shouldn't it be also RGB if nsamples in (3, 4)?

only if image is arranged in Y,X,Band order

[remi-braun]

Transposing before saving correctly saves the colorInterp to RGB (stack = stack.transpose("y", "x", "band")):

Band 1 Block=2048x2048 Type=Byte, ColorInterp=Red
  Description = RED
  NoData Value=0
  Overviews: 10704x6400, 5352x3200, 2676x1600, 1338x800, 669x400
Band 2 Block=2048x2048 Type=Byte, ColorInterp=Green
  Description = GREEN
  NoData Value=0
  Overviews: 10704x6400, 5352x3200, 2676x1600, 1338x800, 669x400
Band 3 Block=2048x2048 Type=Byte, ColorInterp=Blue
  Description = BLUE
  NoData Value=0
  Overviews: 10704x6400, 5352x3200, 2676x1600, 1338x800, 669x400

It could be very valuable to document that the band order should be XYS to work properly with odc.geo.save_cog_with_dask 🙏
(as any pipeline that would use rioxarray or rasterio organize their bands in SYX order)

However, this doesn't solve the nodata issue in the viewer, I think this comes maybe from the missing LAYOUT=COG in the ghost header 🤔

[Kirill888]

@caitlinadams @robbibt ⬆ would be cool to incorporate some of the above to improve docs. I'll be adding documentation label to issues like this and leaving them open for a bit.