Rasterio: WarpedVRT and internal mask

Created on 10 Sep 2018 · 13Comments · Source: mapbox/rasterio

Hello,

First, thank you for this package that makes a lot of gdal power accessible to python developers.

I am encountering an issue I was not able to solve despite reading the multiple issues about similar problems.
I have a cogeotiff created by rasterio. Here's its gdalinfo:

Driver: GTiff/GeoTIFF
Files: 2016-10-29_L8_P185R57.tif
Size is 11376, 11640
Coordinate System is:
PROJCS["WGS 84 / UTM zone 33N",
    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",15],
    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","32633"]]
Origin = (58660.000000000000000,596180.000000000000000)
Pixel Size = (20.000000000000000,-20.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (   58660.000,  596180.000) ( 11d 1'10.28"E,  5d22'50.33"N)
Lower Left  (   58660.000,  363380.000) ( 11d 1'49.82"E,  3d16'46.76"N)
Upper Right (  286180.000,  596180.000) ( 13d 4'13.21"E,  5d23'26.18"N)
Lower Right (  286180.000,  363380.000) ( 13d 4'32.42"E,  3d17' 8.57"N)
Center      (  172420.000,  479780.000) ( 12d 2'56.51"E,  4d20' 5.51"N)
Band 1 Block=256x256 Type=Int16, ColorInterp=Gray
  Overviews: 5688x5820, 2844x2910, 1422x1455, 711x728, 356x364
  Mask Flags: PER_DATASET 
  Overviews of mask band: 5688x5820, 2844x2910, 1422x1455, 711x728, 356x364
Band 2 Block=256x256 Type=Int16, ColorInterp=Undefined
  Overviews: 5688x5820, 2844x2910, 1422x1455, 711x728, 356x364
  Mask Flags: PER_DATASET 
  Overviews of mask band: 5688x5820, 2844x2910, 1422x1455, 711x728, 356x364
Band 3 Block=256x256 Type=Int16, ColorInterp=Undefined
  Overviews: 5688x5820, 2844x2910, 1422x1455, 711x728, 356x364
  Mask Flags: PER_DATASET 
  Overviews of mask band: 5688x5820, 2844x2910, 1422x1455, 711x728, 356x364
Band 4 Block=256x256 Type=Int16, ColorInterp=Undefined
  Overviews: 5688x5820, 2844x2910, 1422x1455, 711x728, 356x364
  Mask Flags: PER_DATASET 
  Overviews of mask band: 5688x5820, 2844x2910, 1422x1455, 711x728, 356x364
Band 5 Block=256x256 Type=Int16, ColorInterp=Undefined
  Overviews: 5688x5820, 2844x2910, 1422x1455, 711x728, 356x364
  Mask Flags: PER_DATASET 
  Overviews of mask band: 5688x5820, 2844x2910, 1422x1455, 711x728, 356x364
Band 6 Block=256x256 Type=Int16, ColorInterp=Undefined
  Overviews: 5688x5820, 2844x2910, 1422x1455, 711x728, 356x364
  Mask Flags: PER_DATASET 
  Overviews of mask band: 5688x5820, 2844x2910, 1422x1455, 711x728, 356x364
Band 7 Block=256x256 Type=Int16, ColorInterp=Undefined
  Overviews: 5688x5820, 2844x2910, 1422x1455, 711x728, 356x364
  Mask Flags: PER_DATASET 
  Overviews of mask band: 5688x5820, 2844x2910, 1422x1455, 711x728, 356x364
Band 8 Block=256x256 Type=Int16, ColorInterp=Undefined
  Overviews: 5688x5820, 2844x2910, 1422x1455, 711x728, 356x364
  Mask Flags: PER_DATASET 
  Overviews of mask band: 5688x5820, 2844x2910, 1422x1455, 711x728, 356x364
Band 9 Block=256x256 Type=Int16, ColorInterp=Undefined
  Overviews: 5688x5820, 2844x2910, 1422x1455, 711x728, 356x364
  Mask Flags: PER_DATASET 
  Overviews of mask band: 5688x5820, 2844x2910, 1422x1455, 711x728, 356x364
Band 10 Block=256x256 Type=Int16, ColorInterp=Undefined
  Overviews: 5688x5820, 2844x2910, 1422x1455, 711x728, 356x364
  Mask Flags: PER_DATASET 
  Overviews of mask band: 5688x5820, 2844x2910, 1422x1455, 711x728, 356x364
Band 11 Block=256x256 Type=Int16, ColorInterp=Undefined
  Overviews: 5688x5820, 2844x2910, 1422x1455, 711x728, 356x364
  Mask Flags: PER_DATASET 
  Overviews of mask band: 5688x5820, 2844x2910, 1422x1455, 711x728, 356x364
Band 12 Block=256x256 Type=Int16, ColorInterp=Undefined
  Overviews: 5688x5820, 2844x2910, 1422x1455, 711x728, 356x364
  Mask Flags: PER_DATASET 
  Overviews of mask band: 5688x5820, 2844x2910, 1422x1455, 711x728, 356x364

It consists in several data bands + an internal mask band.

Steps to reproduce the problem.

I would like to read the internal mask using WarpedVRT.

# profile for vrt
profile = {
    'transform': web_mercator.tile_transform(a=3999, b=4357),
    'height': 256,
    'width': 256,
    'crs': '+init=epsg:3857',
    'resampling': rasterio.enums.Resampling.nearest
}

# open dataset
with rasterio.open(raster) as src:
    print(src.dataset_mask().mean())
    print(src.mask_flag_enums)
    print('---')

    with rasterio.vrt.WarpedVRT(src, **profile) as vrt:
        print(vrt.dataset_mask().mean())
        print(vrt.mask_flag_enums)

Outputs :

77.83712975952267
([<MaskFlags.per_dataset: 2>], [<MaskFlags.per_dataset: 2>], [<MaskFlags.per_dataset: 2>], [<MaskFlags.per_dataset: 2>], [<MaskFlags.per_dataset: 2>], [<MaskFlags.per_dataset: 2>], [<MaskFlags.per_dataset: 2>], [<MaskFlags.per_dataset: 2>], [<MaskFlags.per_dataset: 2>], [<MaskFlags.per_dataset: 2>], [<MaskFlags.per_dataset: 2>], [<MaskFlags.per_dataset: 2>])
---
255.0
([<MaskFlags.all_valid: 1>], [<MaskFlags.all_valid: 1>], [<MaskFlags.all_valid: 1>], [<MaskFlags.all_valid: 1>], [<MaskFlags.all_valid: 1>], [<MaskFlags.all_valid: 1>], [<MaskFlags.all_valid: 1>], [<MaskFlags.all_valid: 1>], [<MaskFlags.all_valid: 1>], [<MaskFlags.all_valid: 1>], [<MaskFlags.all_valid: 1>], [<MaskFlags.all_valid: 1>], [<MaskFlags.all_valid: 1>])

I tried to tweak several parameters (add_alpha, nodata, etc.) but cannot make it work: the vrt mask flag always is all_valid.

Thank you

Versions

Ubuntu 18.04.1
rasterio 1.0.3.post1
gdal 2.2.3

GDAL upstream

Source

thomascoquet

Most helpful comment

yes, GDALRasterBand::GetMaskBand() only returns the alpha band as the mask band in a few particular cases. I'm not sure why this restriction is done. It could probably be removed. In the meantime rasterio could probably also loop through bands and return bands tagged with alpha color interpretation as mask band

rouault on 11 Sep 2018

👍2

All 13 comments

👋 @thomascoquet can you try adding add_alpha=True to profile ?

 with rasterio.open('test.tif') as src:
     print(src.dataset_mask().mean())
     with rasterio.vrt.WarpedVRT(src, add_alpha=True) as dsrc:
         print(dsrc.dataset_mask().mean())

16.243387386202812
16.243387386202812

Edit: oh just realized I'm on GDAL 2.3.1 so it might be different

vincentsarago on 10 Sep 2018

I've already tried that and it produces the same output.
Using your code snipet:

with rasterio.open(raster) as src:
    print(src.dataset_mask().mean())
    print(src.mask_flag_enums)
    print('---')

    with rasterio.vrt.WarpedVRT(src, add_alpha=True) as vrt:
        print(vrt.dataset_mask().mean())
        print(vrt.mask_flag_enums)

outputs:

77.83712975952267
([<MaskFlags.per_dataset: 2>], [<MaskFlags.per_dataset: 2>], [<MaskFlags.per_dataset: 2>], [<MaskFlags.per_dataset: 2>], [<MaskFlags.per_dataset: 2>], [<MaskFlags.per_dataset: 2>], [<MaskFlags.per_dataset: 2>], [<MaskFlags.per_dataset: 2>], [<MaskFlags.per_dataset: 2>], [<MaskFlags.per_dataset: 2>], [<MaskFlags.per_dataset: 2>], [<MaskFlags.per_dataset: 2>])
---
255.0
([<MaskFlags.all_valid: 1>], [<MaskFlags.all_valid: 1>], [<MaskFlags.all_valid: 1>], [<MaskFlags.all_valid: 1>], [<MaskFlags.all_valid: 1>], [<MaskFlags.all_valid: 1>], [<MaskFlags.all_valid: 1>], [<MaskFlags.all_valid: 1>], [<MaskFlags.all_valid: 1>], [<MaskFlags.all_valid: 1>], [<MaskFlags.all_valid: 1>], [<MaskFlags.all_valid: 1>], [<MaskFlags.all_valid: 1>])

thomascoquet on 10 Sep 2018

👍1

@thomascoquet would you be able to share your dataset or a similar file ?

vincentsarago on 10 Sep 2018

@thomascoquet like @vincentsarago I can't reproduce this using GDAL 2.2.4 or 2.3.1. I think it's very possible that GDAL changeset https://trac.osgeo.org/gdal/changeset/41753 has fixed the problem. Can you try upgrading GDAL?

@vincentsarago is also right about adding add_alpha=True – otherwise there is no way for GDAL's VRT code to make any kind of per-dataset mask.

sgillies on 10 Sep 2018

@vincentsarago @sgillies Finally managed to update my GDAL.

Here's my code snippet:

print(rasterio.__version__)
print(rasterio.__gdal_version__)

# open dataset
with rasterio.open(raster) as src:
    print(src.dataset_mask().mean())
    print(src.mask_flag_enums)
    print('---')

    with rasterio.vrt.WarpedVRT(src, add_alpha=True) as vrt:
        print(vrt.dataset_mask().mean())
        print(vrt.mask_flag_enums)

Yields the same result unfortunately:

1.0.3.post1
2.3.1
30.953807830810547
([<MaskFlags.per_dataset: 2>], [<MaskFlags.per_dataset: 2>], [<MaskFlags.per_dataset: 2>], [<MaskFlags.per_dataset: 2>], [<MaskFlags.per_dataset: 2>], [<MaskFlags.per_dataset: 2>], [<MaskFlags.per_dataset: 2>], [<MaskFlags.per_dataset: 2>], [<MaskFlags.per_dataset: 2>], [<MaskFlags.per_dataset: 2>], [<MaskFlags.per_dataset: 2>], [<MaskFlags.per_dataset: 2>])
---
255.0
([<MaskFlags.all_valid: 1>], [<MaskFlags.all_valid: 1>], [<MaskFlags.all_valid: 1>], [<MaskFlags.all_valid: 1>], [<MaskFlags.all_valid: 1>], [<MaskFlags.all_valid: 1>], [<MaskFlags.all_valid: 1>], [<MaskFlags.all_valid: 1>], [<MaskFlags.all_valid: 1>], [<MaskFlags.all_valid: 1>], [<MaskFlags.all_valid: 1>], [<MaskFlags.all_valid: 1>], [<MaskFlags.all_valid: 1>])

Link to small raster generated using gdal_translate + srcwin & geotiff options for cogeo + internal mask.

Here's the gdal info on this file:

Driver: GTiff/GeoTIFF
Files: /media/tom/Johor/update/extract.tif
Size is 512, 512
Coordinate System is:
PROJCS["WGS 84 / UTM zone 33N",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",15],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AUTHORITY["EPSG","32633"]]
Origin = (158660.000000000000000,496180.000000000000000)
Pixel Size = (20.000000000000000,-20.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (  158660.000,  496180.000) ( 11d55'28.50"E,  4d28'57.08"N)
Lower Left  (  158660.000,  485940.000) ( 11d55'29.88"E,  4d23'24.06"N)
Upper Right (  168900.000,  496180.000) ( 12d 1' 0.33"E,  4d28'58.46"N)
Lower Right (  168900.000,  485940.000) ( 12d 1' 1.67"E,  4d23'25.41"N)
Center      (  163780.000,  491060.000) ( 11d58'15.10"E,  4d26'11.26"N)
Band 1 Block=256x256 Type=Int16, ColorInterp=Gray
  Overviews: 256x256, 128x128
  Mask Flags: PER_DATASET 
  Overviews of mask band: 256x256, 128x128
  Metadata:
    name=TRUE_COLOR_R
Band 2 Block=256x256 Type=Int16, ColorInterp=Undefined
  Overviews: 256x256, 128x128
  Mask Flags: PER_DATASET 
  Overviews of mask band: 256x256, 128x128
  Metadata:
    name=TRUE_COLOR_G
Band 3 Block=256x256 Type=Int16, ColorInterp=Undefined
  Overviews: 256x256, 128x128
  Mask Flags: PER_DATASET 
  Overviews of mask band: 256x256, 128x128
  Metadata:
    name=TRUE_COLOR_B
Band 4 Block=256x256 Type=Int16, ColorInterp=Undefined
  Overviews: 256x256, 128x128
  Mask Flags: PER_DATASET 
  Overviews of mask band: 256x256, 128x128
  Metadata:
    name=RDEHZ_G
Band 5 Block=256x256 Type=Int16, ColorInterp=Undefined
  Overviews: 256x256, 128x128
  Mask Flags: PER_DATASET 
  Overviews of mask band: 256x256, 128x128
  Metadata:
    name=RDEHZ_R
Band 6 Block=256x256 Type=Int16, ColorInterp=Undefined
  Overviews: 256x256, 128x128
  Mask Flags: PER_DATASET 
  Overviews of mask band: 256x256, 128x128
  Metadata:
    name=RDEHZ_N
Band 7 Block=256x256 Type=Int16, ColorInterp=Undefined
  Overviews: 256x256, 128x128
  Mask Flags: PER_DATASET 
  Overviews of mask band: 256x256, 128x128
  Metadata:
    name=RDEHZ_S1
Band 8 Block=256x256 Type=Int16, ColorInterp=Undefined
  Overviews: 256x256, 128x128
  Mask Flags: PER_DATASET 
  Overviews of mask band: 256x256, 128x128
  Metadata:
    name=RDEHZ_S2
Band 9 Block=256x256 Type=Int16, ColorInterp=Undefined
  Overviews: 256x256, 128x128
  Mask Flags: PER_DATASET 
  Overviews of mask band: 256x256, 128x128
  Metadata:
    name=NDFI
Band 10 Block=256x256 Type=Int16, ColorInterp=Undefined
  Overviews: 256x256, 128x128
  Mask Flags: PER_DATASET 
  Overviews of mask band: 256x256, 128x128
  Metadata:
    name=LAI
Band 11 Block=256x256 Type=Int16, ColorInterp=Undefined
  Overviews: 256x256, 128x128
  Mask Flags: PER_DATASET 
  Overviews of mask band: 256x256, 128x128
Band 12 Block=256x256 Type=Int16, ColorInterp=Undefined
  Overviews: 256x256, 128x128
  Mask Flags: PER_DATASET 
  Overviews of mask band: 256x256, 128x128

Thank you !

thomascoquet on 11 Sep 2018

with rasterio.open('extract.tif') as src:
    print(src.dataset_mask().mean())
    with rasterio.vrt.WarpedVRT(src, add_alpha=True) as dsrc:
    print(dsrc.dataset_mask().mean())

30.953807830810547
255.0

confirmed too.

thanks @thomascoquet I'm going to deep dive on this.

vincentsarago on 11 Sep 2018

Alright I think I have a good understanding.
This is happening because the dataset has more than 3 bands

# 3 Bands dataset
with rasterio.open('test_rgb.tif') as src:
    print(src.profile)
    print(src.dataset_mask().mean())
    print(src.colorinterp)
    with rasterio.vrt.WarpedVRT(src, add_alpha=True) as dsrc:
        print(dsrc.dataset_mask().mean())

{'driver': 'GTiff', 'dtype': 'uint8', 'nodata': None, 'width': 8192, 'height': 8192, 'count': 3, 'crs': CRS({'init': 'epsg:3857'}), 'transform': Affine(0.5971640348434448, 0.0, -11985326.035115635,
       0.0, -0.5971640348434448, 6574807.4249777235), 'blockxsize': 512, 'blockysize': 512, 'tiled': True, 'compress': 'lzw', 'interleave': 'pixel'}
(<ColorInterp.red: 3>, <ColorInterp.green: 4>, <ColorInterp.blue: 5>)

16.243387386202812
16.243387386202812

# > 3  bands 
with rasterio.open('test.tif') as src:
    print(src.profile)
    print(src.dataset_mask().mean())
    print(src.colorinterp)
    with rasterio.vrt.WarpedVRT(src, add_alpha=True) as dsrc:
        print(dsrc.dataset_mask().mean())

{'driver': 'GTiff', 'dtype': 'uint8', 'nodata': None, 'width': 8192, 'height': 8192, 'count': 9, 'crs': CRS({'init': 'epsg:3857'}), 'transform': Affine(0.5971640348434448, 0.0, -11985326.035115635,
       0.0, -0.5971640348434448, 6574807.4249777235), 'blockxsize': 512, 'blockysize': 512, 'tiled': True, 'compress': 'lzw', 'interleave': 'pixel'}
(<ColorInterp.grey: 1>, <ColorInterp.undefined: 0>, <ColorInterp.undefined: 0>, <ColorInterp.undefined: 0>, <ColorInterp.undefined: 0>, <ColorInterp.undefined: 0>, <ColorInterp.undefined: 0>, <ColorInterp.undefined: 0>, <ColorInterp.undefined: 0>)

16.243387386202812
255.0

Ref: https://github.com/mapbox/rasterio/issues/1403 (this might be linked)

vincentsarago on 11 Sep 2018

also confirmed the alpha band (mask) is only fetched if the dataset is 3 or 1 band

# 2 Bands | INVALID 
with rasterio.open('test_rg.tif') as src:
    print(src.profile)
    print(src.colorinterp)
    print(src.dataset_mask().mean())
    with rasterio.vrt.WarpedVRT(src, add_alpha=True) as dsrc:
        print(dsrc.dataset_mask().mean())

{'driver': 'GTiff', 'dtype': 'uint8', 'nodata': None, 'width': 8192, 'height': 8192, 'count': 2, 'crs': CRS({'init': 'epsg:3857'}), 'transform': Affine(0.5971640348434448, 0.0, -11985326.035115635,
       0.0, -0.5971640348434448, 6574807.4249777235), 'blockxsize': 512, 'blockysize': 512, 'tiled': True, 'compress': 'lzw', 'interleave': 'pixel'}
(<ColorInterp.grey: 1>, <ColorInterp.undefined: 0>)

16.243387386202812
255.0

---


# 1 Band | VALID 
with rasterio.open('test_r.tif') as src:
    print(src.profile)
    print(src.colorinterp)
    print(src.dataset_mask().mean())
    with rasterio.vrt.WarpedVRT(src, add_alpha=True) as dsrc:
        print(dsrc.dataset_mask().mean())

{'driver': 'GTiff', 'dtype': 'uint8', 'nodata': None, 'width': 8192, 'height': 8192, 'count': 1, 'crs': CRS({'init': 'epsg:3857'}), 'transform': Affine(0.5971640348434448, 0.0, -11985326.035115635,
       0.0, -0.5971640348434448, 6574807.4249777235), 'blockxsize': 512, 'blockysize': 512, 'tiled': True, 'compress': 'lzw', 'interleave': 'band'}
(<ColorInterp.grey: 1>,)

16.243387386202812
16.243387386202812

This is a GDAL problem, code that need to be fixed (if possible) is here: https://github.com/OSGeo/gdal/blob/master/gdal/gcore/gdalrasterband.cpp#L6067-L6144

vincentsarago on 11 Sep 2018

@vincentsarago Not sure if it helps, but I have the feeling that if I give src_nodata (without add_alpha), vrt.read(i, masked=True) is correctly masked according to the internal mask (+ eventually data equals to src_nodata).

thomascoquet on 11 Sep 2018

@vincentsarago

with rasterio.open(raster) as src:
    print(src.dataset_mask().mean())
    with rasterio.vrt.WarpedVRT(src, src_nodata=-10001) as vrt:
        print(vrt.dataset_mask().mean())

yields

30.953807830810547
30.953807830810547

Actually, it works only and only if src_nodata value does not exist in raster. If src_nodata = 0, I get:

30.953807830810547
255.0

thomascoquet on 11 Sep 2018

I'm tagging @rouault because he might have a better understanding. To me it seems that GDAL only check for alpha band (when we use add_alpha=True) only for 1 or 3 bands datasets (https://github.com/mapbox/rasterio/issues/1454#issuecomment-420279068).

The src_nodata solution seems to be more a hack that a viable solution IMO.

vincentsarago on 11 Sep 2018

rouault on 11 Sep 2018

👍2

I'm closing this one and marking it an upstream issue. Whether or not we work around as suggested will be tracked in #1456.

sgillies on 11 Sep 2018

Was this page helpful?

0 / 5 - 0 ratings

Related issues

Optimize for "cloud optimized GeoTIFFs"

sgillies · 4Comments

`dataset_mask` return data as int64

vincentsarago · 3Comments

1.1.4 release

sgillies · 3Comments

reproject not taking ndarray indexing into account

davidbrochart · 5Comments

Rasterio 1.0a10

sgillies · 3Comments