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?
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.