1

I know very little of geospatial data. I have a GeoTIFF with a single band that contains a value per pixel (values go from 1 to 6). Using latitude and longitude, I would like to find (i) the coordinates of the pixel and the (ii) value of the relative value.

This is what I could put together so far

dataset  = rio.open(tif)

long = 151.7771477
lat = -32.89674992

row, col = dataset.index(long, lat)
band = dataset.read(1)
band[row][col]

but I get the following error

MemoryError: Unable to allocate 17.4 GiB for an array with shape (1, 125910, 148584) and data type int8
2
  • rasterio seems to try to read the whole geotiff into memory and fails. In this answer to a similar question gis.stackexchange.com/a/368875/15183 window- or block-mode is mentioned. perhaps that helps?
    – til_b
    Mar 16, 2022 at 9:04
  • Thanks - I am looking at it - would you know how to combine the block window option and look for a specific pixel? I am very new and I am not sure how to combine both solutions. Mar 16, 2022 at 9:09

2 Answers 2

4

You need to use the Windowed reading functionality of rasterio as described here

    with rasterio.open(image_file_name) as src:
        meta = src.meta

        # My target coordinates
        x_coord = -110.8270
        y_coord = 32.2215
    
        # Use the transform in the metadata and your coordinates
        rowcol = rasterio.transform.rowcol(meta['transform'], xs=x_coord, ys=y_coord)

        w = src.read(1, window=Window(rowcol[0], rowcol[1], 1, 1))
        # Do something with your data here.
2

You can try working with the transform which can be done using rasterio too. You can put a list of values or just a single number value which I did in the sample code below.

Relevant documentation: https://rasterio.readthedocs.io/en/latest/api/rasterio.transform.html#rasterio.transform.rowcol

Sample code:

import rasterio
from rasterio.windows import Window

image = 'raster.tif'
with rasterio.open(image) as f:
    # Load metadata
    meta = f.meta
    
    # My target coordinates
    x_coord = -110.8270
    y_coord = 32.2215
    
    # Use the transform in the metadata and your coordinates
    rowcol = rasterio.transform.rowcol(meta['transform'], xs=x_coord, ys=y_coord)

    # rowcol value: (977994, 978126)

    y = rowcol[0]
    x = rowcol[1]

    # Load specific pixel only using a window
    window = Window(x,y,1,1)
    arr = f.read(window=window)

print(arr)
# Result is in a nested array
# [[[25281]]]
4
  • this is great but how can I return the value of the pixel? Mar 16, 2022 at 12:17
  • 1
    Sorry, I forgot to include that part of your question. I updated my answer to include the code to get the pixel value. I double checked the pixel values in QGIS and it returns the correct value for me.
    – Esparko
    Mar 16, 2022 at 13:39
  • But will this also read the whole file into memory? I think it will, and thus will run into the same error as in the question. Maybe we can hack something together with windowed reads. I'll give it a try tonight...
    – til_b
    Mar 16, 2022 at 16:57
  • 1
    Hi, sorry for the delay, I updated the code to include using a Window to load a subset of the raster instead of the entire array. Hopefully it works for you.
    – Esparko
    Mar 20, 2022 at 9:17

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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