Rasterio: JPEG 2000 mask issue

[drnextgis]

If read JPEG 2000 image with window that bigger than target image then we get mask=False for output raster. The same operation works as expected for GeoTIFF (returns correct value of mask). Jupyter Notebook with demo: https://gist.github.com/drnextgis/f71038db372946fa81c55836bfc1986f.

Example of output raster:

masked_array(
  data=[[[ 852,  814,  877, ..., 1314, 1408, 1336],
         [ 884,  793,  959, ..., 1344, 1438, 1357],
         [ 876,  878, 1018, ..., 1346, 1404, 1351],
         ...,
         [   0,    0,    0, ...,    0,    0,    0],
         [   0,    0,    0, ...,    0,    0,    0],
         [   0,    0,    0, ...,    0,    0,    0]]],
  mask=False,
  fill_value=999999,
  dtype=uint16)

[drnextgis]

Actually format doesn't matter. It is possible to get the same result for GeoTIFF too. The root cause is here: https://github.com/mapbox/rasterio/blob/fd412a6ea6f21f2d92e06091daad496425f690d8/rasterio/_io.pyx#L393-L394
Does it make sense to change this behavior and masks all pixels that extend beyond the dataset's extent in case of boundless=True and masked=True?

[drnextgis]

The original issue is not fixed. Here is my example:

In [1]: import rasterio

In [2]: from rasterio.windows import Window

In [3]: raster = rasterio.open('/vsicurl/https://transfer.sh/U97tF/rgb.jp2')

In [4]: raster.read(1, boundless=True, masked=True, window=Window(-1, -1, 3, 3))
Out[4]: 
masked_array(
  data=[[0, 0, 0],
        [0, 83, 75],
        [0, 83, 83]],
  mask=[[False, False, False],
        [False, False, False],
        [False, False, False]],
  fill_value=999999,
  dtype=uint8)

[sgillies]

@drnextgis I can't reproduce the problem because I can't access https://transfer.sh/U97tF/rgb.jp2. But locally I've converted the red.tif testing file used for https://github.com/mapbox/rasterio/commit/f627e146d457ad9d6232eff75d49cdc9a0c91508#diff-907f6ae49e082d857f4e4d1f9ce31fbcR64 to a JP2 (using the JP2OpenJPEG driver) and get the mask I'm expecting.

$ rio insp /tmp/red.jp2
Rasterio 1.0.7 Interactive Inspector (Python 3.6.6)
Type "src.meta", "src.read(1)", or "help(src)" for more information.
>>> from rasterio.windows import Window
>>> masked = src.read(1, boundless=True, masked=True, window=Window(-1, -1, 66, 66))
>>> masked
masked_array(
  data=[[--, --, --, ..., --, --, --],
        [--, 204, 204, ..., 204, 204, --],
        [--, 204, 204, ..., 204, 204, --],
        ...,
        [--, 204, 204, ..., 204, 204, --],
        [--, 204, 204, ..., 204, 204, --],
        [--, --, --, ..., --, --, --]],
  mask=[[ True,  True,  True, ...,  True,  True,  True],
        [ True, False, False, ..., False, False,  True],
        [ True, False, False, ..., False, False,  True],
        ...,
        [ True, False, False, ..., False, False,  True],
        [ True, False, False, ..., False, False,  True],
        [ True,  True,  True, ...,  True,  True,  True]],
  fill_value=0,
  dtype=uint8)

My red.jp2 file has a nodata value of 0. What is the case for your file?

[drnextgis]

$ gdalinfo rgb.jp2
Driver: JP2OpenJPEG/JPEG-2000 driver based on OpenJPEG library
Files: rgb.jp2
Size is 100, 100
Coordinate System is:
PROJCS["WGS 84 / Pseudo-Mercator",
    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["Mercator_1SP"],
    PARAMETER["central_meridian",0],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["X",EAST],
    AXIS["Y",NORTH],
    EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs"],
    AUTHORITY["EPSG","3857"]]
Origin = (-6574807.424977720715106,-4070118.882129065692425)
Pixel Size = (76.437028285175984,-76.437028285176893)
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (-6574807.425,-4070118.882) ( 59d 3'45.00"W, 34d18'25.72"S)
Lower Left  (-6574807.425,-4077762.585) ( 59d 3'45.00"W, 34d21'49.84"S)
Upper Right (-6567163.722,-4070118.882) ( 58d59'37.81"W, 34d18'25.72"S)
Lower Right (-6567163.722,-4077762.585) ( 58d59'37.81"W, 34d21'49.84"S)
Center      (-6570985.574,-4073940.734) ( 59d 1'41.40"W, 34d20' 7.79"S)
Band 1 Block=100x100 Type=Byte, ColorInterp=Red
  Overviews: arbitrary
  Image Structure Metadata:
    COMPRESSION=JPEG2000
Band 2 Block=100x100 Type=Byte, ColorInterp=Green
  Overviews: arbitrary
  Image Structure Metadata:
    COMPRESSION=JPEG2000
Band 3 Block=100x100 Type=Byte, ColorInterp=Blue
  Overviews: arbitrary
  Image Structure Metadata:
    COMPRESSION=JPEG2000

rgb.tar.gz

Actually my file doesn't contain nodata value.

[sgillies]

@drnextgis ah, I think I understand now. You'd like this masked read to mask out the area beyond the dataset extents even if the data has no mask and no nodata value?

[drnextgis]

Yes, exactly!

[sgillies]

@drnextgis thanks for your patience! I think #1485 should cover this situation.

[drnextgis]

Thank you Sean! It works now:

In [4]: raster.read(1, boundless=True, masked=True, window=Window(-1, -1, 3, 3))
ERROR 4: No such file or directory
<VRTDataset rasterXSize="101" rasterYSize="101"><SRS>PROJCS["WGS 84 / Pseudo-Mercator",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["Mercator_1SP"],PARAMETER["central_meridian",0],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["X",EAST],AXIS["Y",NORTH],EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs"],AUTHORITY["EPSG","3857"]]</SRS><GeoTransform>-6574883.862006006,76.43702828517598,0.0,-4070042.4451007806,0.0,-76.4370282851769</GeoTransform><VRTRasterBand band="1" dataType="Byte"><ColorInterp>Red</ColorInterp><SimpleSource><SourceFilename relativeToVRT="0">/vsimem/blank-b5a3d9bc-abe7-4762-a65d-b8827a0f6bcd.tif</SourceFilename><SourceBand>1</SourceBand><SourceProperties BlockXSize="3" BlockYSize="3" RasterXSize="101" RasterYSize="101" dataType="Byte" /><SrcRect xOff="0" xSize="3" yOff="0" ySize="3" /><DstRect xOff="0.9999999999992147" xSize="100.0" yOff="0.9999999999992027" ySize="100.00000000000001" /></SimpleSource></VRTRasterBand><VRTRasterBand band="2" dataType="Byte"><ColorInterp>Green</ColorInterp><SimpleSource><SourceFilename relativeToVRT="0">/vsimem/blank-b5a3d9bc-abe7-4762-a65d-b8827a0f6bcd.tif</SourceFilename><SourceBand>2</SourceBand><SourceProperties BlockXSize="3" BlockYSize="3" RasterXSize="101" RasterYSize="101" dataType="Byte" /><SrcRect xOff="0" xSize="3" yOff="0" ySize="3" /><DstRect xOff="0.9999999999992147" xSize="100.0" yOff="0.9999999999992027" ySize="100.00000000000001" /></SimpleSource></VRTRasterBand><VRTRasterBand band="3" dataType="Byte"><ColorInterp>Blue</ColorInterp><SimpleSource><SourceFilename relativeToVRT="0">/vsimem/blank-b5a3d9bc-abe7-4762-a65d-b8827a0f6bcd.tif</SourceFilename><SourceBand>3</SourceBand><SourceProperties BlockXSize="3" BlockYSize="3" RasterXSize="101" RasterYSize="101" dataType="Byte" /><SrcRect xOff="0" xSize="3" yOff="0" ySize="3" /><DstRect xOff="0.9999999999992147" xSize="100.0" yOff="0.9999999999992027" ySize="100.00000000000001" /></SimpleSource></VRTRasterBand></VRTDataset>
Out[4]: 
masked_array(
  data=[[--, --, --],
        [--, 83, 75],
        [--, 83, 83]],
  mask=[[ True,  True,  True],
        [ True, False, False],
        [ True, False, False]],
  fill_value=999999,
  dtype=uint8)

[drnextgis]

What does it mean ERROR 4: No such file or directory?

[drnextgis]

@sgillies is it hard to implement the same behavior in case of boundless warping? I have made a small notebook which illustrates the current state: https://colab.research.google.com/drive/1mPC9GTnHTDDw0Z6x4ti9Bfw1R6_mCpbB

[sgillies]

@drnextgis I think it would be a matter of enhancing the rio-warp command so that you can exploit all the options of warp.reproject(). Would you like to write a ticket for this?