1

I need to search across ~100 1200x1200 raster tiles and produce a min and a max raster. Along the way i also want to exclude values that are over 1 or less than 0. What would be the quickest method to do this in python with gdal? I was previously loading all the rasters into memory and then looping over the pixels, but this seems the slowest possible method.

Thanks!

  • try to use the matrix (extract it from the raster) and look for optimization methods to find the max and min value of big matrix in python. – Gago-Silva Aug 13 '13 at 11:42
1

You could use numpy max/min methods and a masked array.

For example:

from osgeo import gdal
import numpy
import glob

#Loop through and open all rasters and stack them into a 3d array
rasterpaths=glob.glob(r'*.tif')
for ras in rasterpaths:
    ds=gdal.Open(ras)
    dat=ds.GetRasterBand(1).ReadAsArray()
    dat=dat[None,:,:]               #turn 2d array into 3d

    try:stk=numpy.vstack((stk,dat)) #Do we already have a 3d stack?
    except NameError:stk=dat        #Nope, this is the first time through the loop

#Create masked array where valid values are >=0 and <=1
stk=numpy.ma.masked_outside(stk,0,1)

#Get max and min
max=stk.max(axis=0)
min=stk.min(axis=0)

#Write out to new rasters with gdal if you like...
0

For anyone in the future: I solved this by essentially only loading two bands into memory (i and i+1) and then comparing just these two. If i+1 met was greater/smaller it was accepted as the new maximum/minimum.

The theory is in very similar to https://stackoverflow.com/questions/11964450/python-order-a-list-of-numbers-without-built-in-sort-min-max-function

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.