Different values between a windowed read and a normal read in rasterio

Question

I am struggling to understand why I obtain different results using the same coordinates in this example:

import rasterio
#data link: https://earthexplorer.usgs.gov/download/external/options/SENTINEL_2A/3022286/INVSVC/

pathToImgFolder = ""
pathData = pathToImgFolder + "L1C_T21LXC_A001666_20170701T140052\S2B_MSIL1C_20170701T140049_N0205_R067_T21LXC_20170701T140052.SAFE\GRANULE\L1C_T21LXC_A001666_20170701T140052\IMG_DATA\T21LXC_20170701T140049_B01.jp2" 
boudingBoxCoordinates = (606859.0750363453, 8241169.269917269, 607219.0750363453, 8241529.269917269)
points = [(607014.0750363453,8241374.26991727),(607024.0750363444,8241374.26991727),\
          (607034.0750363453,8241374.26991727),(607044.0750363453,8241374.26991727),\
          (607054.0750363442,8241374.26991727)]

with rasterio.open(pathData) as img:
    bandWindow = rasterio.windows.from_bounds(*boudingBoxCoordinates, img.transform)
    winTransform = rasterio.windows.transform(bandWindow,img.transform)
    bandData = img.read(1, window = bandWindow)
    allData = img.read(1)
    for point in points:
        rBand,cBand = rasterio.transform.rowcol(winTransform,point[0],point[1])
        rFull,cFull = img.index(point[0],point[1])
        bandVal = bandData[rBand,cBand]
        fullVal = allData[rFull,cFull]
        print(f"{fullVal} {'=' if bandVal == fullVal else '!='} {bandVal}")

When I execute this example i get this result:

1202 = 1202
1330 != 1202
1330 != 1202
1330 = 1330
1330 = 1330

I think it is related to the spatial resolution of the band(60m), but I am not sure and I do not know what to do.

The img.profile is:

{'driver': 'JP2OpenJPEG', 'dtype': 'uint16', 'nodata': None, 'width': 1830, 'height': 1830, 'count': 1, 'crs': CRS.from_epsg(32721), 'transform': Affine(60.0, 0.0, 600000.0, 0.0, -60.0, 8300020.0), 'blockxsize': 192, 'blockysize': 192, 'tiled': True} 

Answer 1

Your window has floating point offsets.

If I use a window with float offsets to read from an image and then write out to a cropped image, the output data is shifted slightly, enough to make a difference. I don't know if this is intentional or a bug.

You will get the expected results if you round your window offsets, using the window.round_offsets() method.

with rasterio.open(pathData) as img:
    bandWindow = rasterio.windows.from_bounds(
        *boudingBoxCoordinates, 
        img.transform).round_offsets()
    etc...

Here's an example of just how much the shift is with the floating point window:

import numpy as np
import rasterio as rio
from rasterio.transform import Affine
from rasterio.crs import CRS

print(rio.__version__)    
# '1.1.5'

# Make some dummy data
profile = {'driver': 'GTIFF', 'dtype': 'uint16', 'nodata': None,
           'width': 1830, 'height': 1830, 'count': 1,
           'crs': CRS.from_epsg(32721), 'transform': Affine(60.0, 0.0, 600000.0, 0.0, -60.0, 8300020.0),
           'blockxsize': 192, 'blockysize': 192, 'tiled': True}

zmin, zmax = 1, 1234
raster = np.random.randint(zmin,zmax, (profile['count'], profile['height'], profile['width'])).astype(profile['dtype'])

with rio.open('/tmp/test.tif', 'w', **profile) as ds:
    ds.write(raster)


bbox = (606859.0750363453, 8241169.269917269, 607219.0750363453, 8241529.269917269)

with rio.open('/tmp/test.tif') as src:
    
    # floating point window
    fwindow = rio.windows.from_bounds(*bbox, src.transform)
    ftransform = src.window_transform(fwindow)

    # rounded window
    iwindow = fwindow.round_offsets()
    itransform = src.window_transform(iwindow)

    # Write out the data with the floating point window
    kwargs = src.meta.copy()
    kwargs.update({
        'height': fwindow.height,
        'width': fwindow.width,
        'transform': ftransform})

    with rio.open('/tmp/croppedf.tif', 'w', **kwargs) as dst:
        dst.write(src.read(window=fwindow))

    # Write out the data with the rounded window
    kwargs = src.meta.copy()
    kwargs.update({
        'height': iwindow.height,
        'width': iwindow.width,
        'transform': itransform})

    with rio.open('/tmp/croppedi.tif', 'w', **kwargs) as dst:
        dst.write(src.read(window=iwindow))