A certain combination of parameters leads to unexpected results

[drnextgis]

I would expect to get the same set of mask_flags for Case 4 as for the rest of the cases. Am I missing something or it is a bug?

$ rio info original.tif | jq -c .mask_flags,.nodata
[["nodata"],["nodata"],["nodata"]]
0

🟢 Case 1

--add-mask

$ rio cogeo create --add-mask original.tif cogeo1.tif
$ rio info cogeo1.tif | jq -c .mask_flags
[["per_dataset"],["per_dataset"],["per_dataset"]]

🟢 Case 2

--add-mask --web-optimized

$ rio cogeo create --add-mask --web-optimized original.tif cogeo2.tif
$ rio info cogeo2.tif | jq -c .mask_flags                           
[["per_dataset"],["per_dataset"],["per_dataset"]]

🟢 Case 3

--add-mask --use-cog-driver

$ rio cogeo create --add-mask --use-cog-driver original.tif cogeo3.tif
$ rio info cogeo3.tif | jq -c .mask_flags  
[["per_dataset"],["per_dataset"],["per_dataset"]]

🔴 Case 4

--add-mask --use-cog-driver --web-optimized

$ rio cogeo create --add-mask --use-cog-driver --web-optimized original.tif cogeo4.tif
$ rio info cogeo4.tif | jq -c .mask_flags 
[["per_dataset","alpha"],["per_dataset","alpha"],["per_dataset","alpha"],["all_valid"]]

[vincentsarago]

Thanks for the report @drnextgis
That's an interesting issue. Going thought the code I can't find anything suspicious. Can you try doing the same with GDAL commands (use gdalwarp -of VRT -dstalpha in.tif out.vrt + gdal_translate -of COG -mask 4 out.vrt out.tif, or something like this)

[drnextgis]

It works fine. At least it doesn't produce an alpha band which is my main problem.

$ gdalinfo --version
GDAL 3.4.0dev-8dddba9d9065bf0a2f41be1282ff2dba79c21797, released 2021/06/16

$ gdalwarp -of VRT -dstalpha original.tif out.vrt
Creating output file that is 8366P x 10042L.
Processing original.tif [1/1] : 0Using internal nodata values (e.g. 0) for image original.tif.
...10...20...30...40...50...60...70...80...90...100 - done.

$ gdal_translate -of COG -mask 4 out.vrt out.tif
Input file size is 8366, 10042
0...10...20...30...40...50...60...70...80...90...100 - done.

$ rio info out.tif | jq -c .mask_flags  
[["per_dataset"],["per_dataset"],["per_dataset"],["per_dataset"]]

[vincentsarago]

can you print the compression used for all the COG ?

[drnextgis]

$ rio info original.tif| jq -c .compress 
null
$ rio info cogeo1.tif | jq -c .compress   
"deflate"
$ rio info cogeo2.tif | jq -c .compress
"deflate"
$ rio info cogeo3.tif | jq -c .compress
"deflate"
$ rio info cogeo4.tif | jq -c .compress
"deflate"
$ rio info out.tif | jq -c .compress
null

[vincentsarago]

@drnextgis is your file in byte ?

[drnextgis]

$ gdalinfo original.tif
Driver: GTiff/GeoTIFF
Files: original.tif
Size is 8366, 10042
Coordinate System is:
PROJCS["WGS 84 / UTM zone 38N",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",45],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["Easting",EAST],
    AXIS["Northing",NORTH],
    AUTHORITY["EPSG","32638"]]
Origin = (615076.000000000000000,3441809.000000000000000)
Pixel Size = (1.000000000000000,-1.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (  615076.000, 3441809.000) ( 46d12'24.08"E, 31d 6'16.19"N)
Lower Left  (  615076.000, 3431767.000) ( 46d12'19.96"E, 31d 0'50.06"N)
Upper Right (  623442.000, 3441809.000) ( 46d17'39.83"E, 31d 6'13.13"N)
Lower Right (  623442.000, 3431767.000) ( 46d17'35.41"E, 31d 0'47.01"N)
Center      (  619259.000, 3436788.000) ( 46d14'59.82"E, 31d 3'31.63"N)
Band 1 Block=8366x1 Type=Byte, ColorInterp=Red
  NoData Value=0
Band 2 Block=8366x1 Type=Byte, ColorInterp=Green
  NoData Value=0
Band 3 Block=8366x1 Type=Byte, ColorInterp=Blue
  NoData Value=0

[drnextgis]

@vincentsarago you can reproduce it with publicly available dataset:

$ rio cogeo create --add-mask --use-cog-driver --web-optimized \
/vsicurl/https://github.com/mapbox/rasterio/blob/master/tests/data/RGB.byte.tif\?raw\=true rasterio.tif
$ rio info rasterio.tif | jq -c .mask_flags 
[["per_dataset","alpha"],["per_dataset","alpha"],["per_dataset","alpha"],["all_valid"]]

[vincentsarago]

🙏 @drnextgis

So this is a GDAL bug

$ gdalwarp -of VRT -dstalpha RGB.byte.tif RGB.byte.vrt 

$ gdal_translate -of COG -mask 4 -co TILING_SCHEME=GoogleMapsCompatible RGB.byte.vrt RGB.cog.tif

$ rio info RGB.cog.tif | jq '.mask_flags' -c 
>>> [["per_dataset","alpha"],["per_dataset","alpha"],["per_dataset","alpha"],["all_valid"]]

The problem is that GDAL will use an Alpha band when retrojecting the data to WebMercator

See https://gdal.org/drivers/raster/cog.html#reprojection-related-creation-options

ADD_ALPHA=YES/NO: Whether an alpha band is added in case of reprojection. Defaults to YES.

So by default, our mask band will be translated to an Alpha Band.

The only way you can keep the mask is to use the JPEG compression (which anyway works better with internal mask than deflate dataset)

$ rio cogeo create --add-mask --web-optimized --use-cog-driver RGB.byte.tif rgb_cog_jpeg.tif -p jpeg            
$ rio info rgb_cog_jpeg.tif | jq '.mask_flags' -c 
[["per_dataset"],["per_dataset"],["per_dataset"]]

What's next

We could add warning if user use --add-mask for gdal cog driver/web-optimized and non JPEG compression here https://github.com/cogeotiff/rio-cogeo/blob/master/rio_cogeo/cogeo.py#L356 🤷‍♂️

[drnextgis]

Thank you @vincentsarago for sorting this out!