Point out of image extent (negative row,col) returned value

[Olyna]

Expected behavior and actual behavior.

Compare 3 different implementations to read pixel value of a raster, using point coordinates. Reading implementations are:

1. use point coordinates with rasterio sample()
2. find image coordinates and then use rasterio windowed read()
3. find image coordinates and then use gdal ReadAsArray()

Both points and raster should have the same Coordinate Reference System. Points do not intersect with raster's tile index. The expected behavior for all points and all implementations is to return None or empty array. The returned pixel value seems to be inconsistent when image coordinates (row, column) are negative.

Implementation 1:
for all points -also positive image coords- pixelValue = [0] (why?)
Implementation 2:
if column number == -1 then pixelValue = empty array ✓
if column number == -2 then pixelValue = raster's no-data value (why?)
if row number == -1 then pixelValue = empty array ✓
if row number == -2 then pixelValue = raster's no-data value (why?)
Implementation 3:
if column number == -1 then pixelValue = raster's no-data value (why?)
if column number == -2 then pixelValue = None ✓
if row number == -1 then pixelValue = raster's no-data value (why?)
if row number == -2 then pixelValue = None ✓

Steps to reproduce the problem.

Data
Download raster image here.
Download points shapefile here.

Implemented functions

import math
import pyproj
import geopandas as gpd
from osgeo import gdal, osr
import rasterio as rio
from rasterio.windows import Window

def rio_point_sample(tif_fp:str, x:float, y:float, xy_crs:pyproj.crs.crs.CRS):

    # Image value using method sample()
    with rio.open(tif_fp, 'r') as src:
        # Image and point CRS must be the same
        assert src.crs == xy_crs
        try:
            cellValue = next(src.sample([(x, y)]))
        except Exception as e:
            cellValue = e
    return cellValue


def rio_point_readwindow(tif_fp:str, x:float, y:float, xy_crs:pyproj.crs.crs.CRS):

    # Image value using window and method read()
    with rio.open(tif_fp, 'r') as src:
        # Image and point CRS must be the same
        assert src.crs == xy_crs
        # Image coordinates
        rowNum, colNum = src.index(x, y)
        # One pixel window
        window = Window.from_slices(
            rows=(rowNum, rowNum+1),
            cols=(colNum, colNum+1),
            height=src.height,
            width=src.width)
        # Read only the first band using window
        try:
            cellValue = src.read(1, window=window)
        except Exception as e:
            cellValue = e
    return cellValue, rowNum, colNum


def gdal_point_readasarray(tif_fp:str, x:float, y:float, xy_crs:pyproj.crs.crs.CRS):

    dataset = gdal.Open(tif_fp)
    # Image and point CRS must be the same
    srs = osr.SpatialReference()
    srs.ImportFromWkt(dataset.GetProjection())
    assert srs.GetAuthorityCode(None) == xy_crs
    # Image coordinates
    transform = dataset.GetGeoTransform()
    colNum = int(math.floor((x - transform[0]) / transform[1]))
    rowNum = int(math.floor((transform[3] - y) / -transform[5]))
    # Read the first band
    band = dataset.GetRasterBand(1)
    # Cell value
    try:
        cellValue = band.ReadAsArray(colNum, rowNum, 1, 1)
    except Exception as e:
        cellValue = e
    return cellValue, rowNum, colNum

Environment Information

$ rio --show-versions
rasterio info:
  rasterio: 1.3.8
      GDAL: 3.7.0
      PROJ: 9.2.1
      GEOS: 3.11.2
 PROJ DATA: ~/anaconda3/envs/pygdal/share/proj
 GDAL DATA: ~/anaconda3/envs/pygdal/share/gdal

System:
    python: 3.11.4 | packaged by conda-forge | (main, Jun 10 2023, 18:08:17) [GCC 12.2.0]
executable: ~/anaconda3/envs/pygdal/bin/python
   machine: Linux-5.19.0-46-generic-x86_64-with-glibc2.35

Python deps:
    affine: 2.4.0
     attrs: 23.1.0
   certifi: 2023.05.07
     click: 8.1.3
     cligj: 0.7.2
    cython: None
     numpy: 1.25.0
    snuggs: 1.4.7
click-plugins: None
setuptools: 68.0.0

Installation Method

conda 4.10.3

[groutr]

Regarding the usage of rasterio.sample():

https://github.com/rasterio/rasterio/blob/main/rasterio/sample.py#

If the coordinate is out of bounds, then the dataset's nodata or 0 is used as the return value. That could explain why you are getting 0. What is the dataset's nodata value?

I'll have to look deeper to understand what is going on with the other rasterio implementation.

[sgillies]

@Olyna thank you for bringing this question!

Rasterio's read() is intended to behave like numpy array slicing. Think of a window as a pair of slice objects. When the window is outside the raster indexes, rasterio returns an empty array. Now, if you do read(..., boundless=True) you will get something different: an array containing nodata values.

The sample() method acts as though you have passed boundless=True. Thus, you should get nodata values for points outside the raster.

If I've read your description correctly, what you observe is to be expected. sample() is consistent with read(..., boundless=True).