[GIS] Reading subset of GeoTIFF file without having to read whole file first

gdalgeotiff-tiffpythonremote sensing

I have a set of GeoTIFF files which I have generated from the MODIS LST Product and have used the coordinate information present within the file to work out the index of my 'pixel of interest'/study site. My aim is to extract the data value at this index and so far have been using this code within python:

myfile = band.ReadAsArray(0, 0, cols, rows)
value = myfile[y_index, x_index]

which works just fine. This generates the result I want but I am having to read in the whole data file first, before obtaining the result I require. Is it possible to read just my 'pixel of interest' or a subset of the data without having to read the whole file into memory beforehand?

I am using Python and the gdal module within it.

Best Answer

I was trying to make things more complicated than they needed to be, one simple solution is as follows:

When using ReadAsArray(), the function works as follows:

data = myfile.ReadAsArray(x_offset, y_offset, x_size, y_size)

So if you use:

data = myfile.ReadAsArray(x_index, y_index, 1, 1)

This will read in JUST one pixel.

But the code can be adapted to read in a 3x3 window etc. using:

data = myfile.ReadAsArray(x_index - 1, y_index - 1, 3, 3)

Related Solutions

Python – Efficiently Reading VRT GeoTIFF Image Stacks

From your question, you appear to want to read time series of individual pixels. The fastest way I have found for this is to convert the VRT to a geotiff with option "INTERLEAVE=BAND".

If that's not an option (because it takes loads of disk space), you can also use a ThreadPool:

from concurrent.futures import ThreadPoolExecutor
from functools import partial
from osgeo import gdal

def read_pxl(band,x,y, g):
    return g.GetRasterBand(band).ReadAsArray(xoff=x, yoff=y,
                                win_xsize=1, win_ysize=1)
g = gdal.Open(my_super_duper_fname)
n_bands = g.RasterCount
wrapper=partial(read_pxl, g=g, x=1800, y=600)
with ThreadPoolExecutor(max_workers=20) as exec:
    retval = exec.map(wrapper, range(1, n_bands + 1))

In your case, if you use the threadpool it may just be faster to index the independent GeoTIFF files directly without going via the VRT.

Best Answer

Related Solutions

Python – Efficiently Reading VRT GeoTIFF Image Stacks

Related Question