1

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} 
  • Hello @user2856, every once in a while the website just did not work, so my 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} – AdrianoII Sep 11 at 3:49
5

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))

enter image description here

| improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.