-1

I have several large rasters (300MB - 1GB), each raster has been reclassified with numpy array and GDAL blocks. I now need to create an output file with only the maximum value per pixel used in the output. If I just try to open each image as an array, I run out of memory. (in ESRI i would use cell statistics set to MAX), but this needs to be accomplished with open source programming.

I am trying to open each file by block:

for i in range(0, rows, yBlocksize):
    if i + yBlocksize < rows:
        numRows = yBlocksize
    else:
        numRows = rows - i

        # loop through columns
    for j in range(0, cols, xBlocksize):
        if j + xBlocksize < cols:
            numCols = xBlocksize
        else:
            numCols = cols - j
        array1 = band1.ReadAsArray(j, i, numCols, numRows)
        print ("array1")
        print array1.shape

        array2 = band2.ReadAsArray(j, i, numCols, numRows)
        print ("array2")
        print array2.shape

        array3 = band3.ReadAsArray(j, i, numCols, numRows)
        print ("array3")
        print array3.shape

        array4 = band4.ReadAsArray(j, i, numCols, numRows)
        print ("array4")
        print array4.shape

From there, I have found several options for finding the 'max' however, I am getting stuck in the 'for' loop and my script will continuously print 'array1, array2' ect. If I add 'break' after array 4 it is ignored.

Other times I get a zero output, or errors: "raise ValueError("expected array of dim 2") ValueError: expected array of dim 2"

I have to be missing something simple. running: python 2.7 GDAL 2.1.3 numpy 1.13.0

I have looked at the following links and tried pretty much all of these options: How to operate maximum value compose of some raster image?

gdal/python WriteArray working for 8-bit but not 16-bit

Optimising Min/Max temporal raster search in python/gdal

create a maximum raster using gdal_calc

I have also looked here - http://tretherington.blogspot.com/2015/04/some-stuff-import-numpy-as-npa-np.html

so it has to be the size of my raster data and/or reading the data into blocks.

  • Shouldn't the calculation for xMax and YMax be independent of each other? two separate for loops, rather than nested – P.T. Curran Aug 1 '17 at 19:30
  • The way the 'for' loop is written is what I have learned through gis.usu.edu/~chrisg/python/2009 and have had great success in reclassifying the image with that format. the 'array' after that, I am not sure about. – als Aug 1 '17 at 19:53
  • You're right, @als. I missed the use of the variable i inside the second for loop. – P.T. Curran Aug 2 '17 at 13:07
  • SOLVED: through a combination of already sited locations, and no function was needed. This solves with reading each array into blocks, – als Aug 7 '17 at 18:56
2

SOLVED: through a combination of already sited locations, and no function was needed. This solves with reading each array into blocks,

for i in   range(0, rows, yBlocksize):
    if i + yBlocksize < rows:
        numRows = yBlocksize
    else:
        numRows = rows - i

        for j in range(0, cols, xBlocksize):
            if j + xBlocksize < cols:
            numCols = xBlocksize
        else:
            numCols = cols - j

        array1 = band1.ReadAsArray(j, i, numCols, numRows)
        array2 = band2.ReadAsArray(j, i, numCols, numRows)
        array3 = band3.ReadAsArray(j, i, numCols, numRows)
        array4 = band4.ReadAsArray(j, i, numCols, numRows)

        z = np.dstack((array1, array2, array3, array4))
        print z.shape
        maxIndex = np.argmax(z, axis=2)
        nRow, nCol = np.shape(array1)
        col, row = np.meshgrid(range(nCol), range(nRow))

        maxValue = z[row, col, maxIndex]
        print maxValue

        outBand2a.WriteArray(maxValue, j, i)

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.