Support for COGs with nodata masks

[gadomski]

I'm trying to render this COG: https://maxar-opendata.s3.amazonaws.com/events/yellowstone-flooding22/ard/12/120000020112/2022-06-18/10300100D51B8C00-visual.tif. I edited COG_URL in the cog-basic example, loaded it up, and got this error:

Error: Unexpected number of channels in raster data: 1
    addAlphaChannel geotiff.ts:103
    getTileData render-pipeline.ts:143

The COG appears to be good, three-band data: https://cogeo.org/map/#/url/https%3A%2F%2Fmaxar-opendata.s3.amazonaws.com%2Fevents%2Fyellowstone-flooding22%2Fard%2F12%2F120000020112%2F2022-06-18%2F10300100D51B8C00-visual.tif/center/-110.5844,45.4453/zoom/14

The `gdal info`

Driver: GTiff/GeoTIFF
Files: /Users/gadomski/Downloads/10300100D51B8C00-visual.tif
Size is 17408, 17408
Coordinate System is:
PROJCRS["WGS 84 / UTM zone 12N",
    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 12N",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",-111,
            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["Navigation and medium accuracy spatial referencing."],
        AREA["Between 114°W and 108°W, northern hemisphere between equator and 84°N, onshore and offshore. Canada - Alberta; Northwest Territories (NWT); Nunavut; Saskatchewan. Mexico. United States (USA)."],
        BBOX[0,-114,84,-108]],
    ID["EPSG",32612]]
Data axis to CRS axis mapping: 1,2
Origin = (529843.750000000000000,5065156.250000000000000)
Pixel Size = (0.305175781250000,-0.305175781250000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  LAYOUT=COG
  SOURCE_COLOR_SPACE=YCbCr
  COMPRESSION=YCbCr JPEG
  INTERLEAVE=PIXEL
  JPEG_QUALITY=75
  JPEGTABLESMODE=1
Corner Coordinates:
Upper Left  (  529843.750, 5065156.250) (110d36'58.97"W, 45d44'21.55"N)
Lower Left  (  529843.750, 5059843.750) (110d37' 0.15"W, 45d41'29.41"N)
Upper Right (  535156.250, 5065156.250) (110d32'53.14"W, 45d44'20.65"N)
Lower Right (  535156.250, 5059843.750) (110d32'54.52"W, 45d41'28.51"N)
Center      (  532500.000, 5062500.000) (110d34'56.69"W, 45d42'55.05"N)
Band 1 Block=512x512 Type=Byte, ColorInterp=Red
  Overviews: 8704x8704, 4352x4352, 2176x2176, 1088x1088, 544x544, 272x272
  Mask Flags: PER_DATASET 
  Overviews of mask band: 8704x8704, 4352x4352, 2176x2176, 1088x1088, 544x544, 272x272
Band 2 Block=512x512 Type=Byte, ColorInterp=Green
  Overviews: 8704x8704, 4352x4352, 2176x2176, 1088x1088, 544x544, 272x272
  Mask Flags: PER_DATASET 
  Overviews of mask band: 8704x8704, 4352x4352, 2176x2176, 1088x1088, 544x544, 272x272
Band 3 Block=512x512 Type=Byte, ColorInterp=Blue
  Overviews: 8704x8704, 4352x4352, 2176x2176, 1088x1088, 544x544, 272x272
  Mask Flags: PER_DATASET 
  Overviews of mask band: 8704x8704, 4352x4352, 2176x2176, 1088x1088, 544x544, 272x272

I tried to do some debugging on my side, but I quickly reached the extent of my knowledge/ability.

[kylebarron]

What's happening is that the geotiff.js readRasters command is only loading one channel. I'm not sure why

[kylebarron]

But then when I try to force reading 3 samples, I get an error

RangeError: Invalid sample index '1'.
    at GeoTIFFImage.readRasters (geotiffimage.js:588:33)
    at getTileData (render-pipeline.ts:131:67)
    at COGLayer._getTileData (cog-layer.ts:307:24)
    at getTileData (cog-layer.ts:428:41)
    at TileLayer.getTileData (tile-layer.ts:327:14)
    at Tile2DHeader._loadData (tile-2d-header.ts:137:24)

[kylebarron]

Here's the rio cogeo info:

> rio cogeo info https://maxar-opendata.s3.amazonaws.com/events/yellowstone-flooding22/ard/12/120000020112/2022-06-18/10300100D51B8C00-visual.tif
Driver: GTiff
File: https://maxar-opendata.s3.amazonaws.com/events/yellowstone-flooding22/ard/12/120000020112/2022-06-18/10300100D51B8C00-visual.tif
COG: True
Compression: YCbCr JPEG
ColorSpace: YCbCr

Profile
    Width:            17408
    Height:           17408
    Bands:            3
    Tiled:            True
    Dtype:            uint8
    NoData:           None
    Alpha Band:       False
    Internal Mask:    True
    Interleave:       PIXEL
    ColorMap:         False
    ColorInterp:      ('red', 'green', 'blue')
    Scales:           (1.0, 1.0, 1.0)
    Offsets:          (0.0, 0.0, 0.0)

Geo
    Crs:              EPSG:32612
    Origin:           (529843.75, 5035156.25)
    Resolution:       (0.30517578125, -0.30517578125)
    BoundingBox:      (529843.75, 5029843.75, 535156.25, 5035156.25)
    MinZoom:          11
    MaxZoom:          18

Image Metadata
    AREA_OR_POINT: Area

Image Structure
    LAYOUT: COG
    SOURCE_COLOR_SPACE: YCbCr
    COMPRESSION: YCbCr JPEG
    INTERLEAVE: PIXEL
    JPEG_QUALITY: 75
    JPEGTABLESMODE: 1

Band 1
    ColorInterp: red

Band 2
    ColorInterp: green

Band 3
    ColorInterp: blue

IFD
    Id      Size           BlockSize     Decimation
    0       17408x17408    512x512       0
    1       8704x8704      512x512       2
    2       4352x4352      512x512       4
    3       2176x2176      512x512       8
    4       1088x1088      512x512       16
    5       544x544        512x512       32
    6       272x272        512x512       64

One thing to note is that in JS the image (IFD) count is 14, which is double what rio cogeo info says:

So I'm guessing that maybe there are nodata masks in this image? I haven't come across a file with nodata masks before, so if that's the case it would make sense

[kylebarron]

Yeah, that looks likely. Here's fetching the data from the first and second IFD. The second IFD looks to be a boolean array:

[kylebarron]

Oh also from the rio cogeo info 😅:

Profile
    Internal Mask:    True

[kylebarron]

Looks like we need to:

• parse all the IFDs
• look for NewSubFileType being not an overview
• PhotometricInterpretation should be 4 in the cast of masks

Ideally we'll start adding tests here from https://github.com/developmentseed/geotiff-test-data for inferring default rendering behavior for different sorts of GeoTIFFs

[gadomski]

developmentseed/geotiff-test-data#25

Not sure if this is fully duplicative on the maxar example?

[kylebarron]

I think #217 will be the way to go towards solving this issue

[kylebarron]

Closing as duplicate of #168