I'm following this guide to extract pixels with certain rgb code (eg 255,0,0) in a .tif raster using Anaconda 32bit:
RasterToNumPyArray—ArcPy Functions | ArcGIS for Desktop
But I receive this error:
os.chdir("C:\\Users\\ivn\\Desktop\\NumpArr") filein = os.path.join(os.getcwd(),r'AR00199.tif') fGDB = "fGDB.gdb" if not os.path.exists(os.path.join(os.getcwd(),fGDB)): arcpy.CreateFileGDB_management(os.getcwd(),fGDB) fileoutgdb = os.path.join(os.getcwd(),r"fGDB.gdb\blocco") arcpy.env.overwriteOutput = True arcpy.env.outputCoordinateSystem = filein arcpy.env.cellSize = filein arcpy.env.workspace = os.path.join(os.getcwd(),r"fGDB.gdb") arcpy.env.scratchWorkspace = os.path.join(os.getcwd(),r"fGDB.gdb") inRas = arcpy.Raster(filein) b1 = arcpy.sa.Raster(r'AR00199.tif\Band_1') b2 = arcpy.sa.Raster(r'AR00199.tif\Band_2') b3 = arcpy.sa.Raster(r'AR00199.tif\Band_3') blocksize = 512 blockno = 0 for x in range(0,inRas.width,blocksize): for y in range(0, inRas.height, blocksize): blockno += 1 mx = inRas.extent.XMin + x * inRas.meanCellWidth my = inRas.extent.YMin + y * inRas.meanCellHeight lx = min([x + blocksize, inRas.width]) ly = min([y + blocksize, inRas.height]) b1nr = arcpy.RasterToNumPyArray(b1,arcpy.Point(mx, my),lx-x, ly-y) b2nr = arcpy.RasterToNumPyArray(b2,arcpy.Point(mx, my),lx-x, ly-y) b3nr = arcpy.RasterToNumPyArray(b3,arcpy.Point(mx, my),lx-x, ly-y) rgb = np.array([b1nr, b2nr,b3nr], np.uint8) #Operations on blocks.... .... .... .... myRasterBlock = arcpy.NumPyArrayToRaster(output_array,\ arcpy.Point(mx, my),inRas.meanCellWidth,inRas.meanCellHeight) filetemp = '_'.join([fileoutgdb,str(blockno)]) myRasterBlock.save(filetemp)
It save about 173 grids with blocksize = 512 and 73 whit blocksize = 1024 in gdb and then I receive the error... why?? Can someone help?
This is the result of the operation for a block:
Solved! Go to Solution.
My guess is something is happening along the way with "#Operations on blocks" and arcpy.NumPyArrayToRaster() may not be returning a valid array. I would do some debugging of myRasterBlock before saving it, e.g., insert some print statements to print out properties of the raster. Also, print out filetemp . If you could include the results of those print statements, it might help people provide more feedback.
My guess is something is happening along the way with "#Operations on blocks" and arcpy.NumPyArrayToRaster() may not be returning a valid array. I would do some debugging of myRasterBlock before saving it, e.g., insert some print statements to print out properties of the raster. Also, print out filetemp . If you could include the results of those print statements, it might help people provide more feedback.
You use the term 'block' ... do you mean memory size block or a physical block as in the number of rows and columns.
Here is a 3D array and its construction
import numpy as np a = np.ones((3,6,6)) # make a 3D array a[0] = np.arange(36).reshape(6,6) # 36 sequential numbers a[1] = 2 # assign all 2 a[2] = 3 # ditto but use 3 >>> a # like an image array([[[ 0., 1., 2., 3., 4., 5.], [ 6., 7., 8., 9., 10., 11.], [ 12., 13., 14., 15., 16., 17.], [ 18., 19., 20., 21., 22., 23.], [ 24., 25., 26., 27., 28., 29.], [ 30., 31., 32., 33., 34., 35.]], [[ 2., 2., 2., 2., 2., 2.], [ 2., 2., 2., 2., 2., 2.], [ 2., 2., 2., 2., 2., 2.], [ 2., 2., 2., 2., 2., 2.], [ 2., 2., 2., 2., 2., 2.], [ 2., 2., 2., 2., 2., 2.]], [[ 3., 3., 3., 3., 3., 3.], [ 3., 3., 3., 3., 3., 3.], [ 3., 3., 3., 3., 3., 3.], [ 3., 3., 3., 3., 3., 3.], [ 3., 3., 3., 3., 3., 3.], [ 3., 3., 3., 3., 3., 3.]]])
Now you can block it in 3 or 2 dimensions etc. This is the top left 3x3 block for all 3 'bands'
>>> a[:,:3,:3] array([[[ 0., 1., 2.], [ 6., 7., 8.], [ 12., 13., 14.]], [[ 2., 2., 2.], [ 2., 2., 2.], [ 2., 2., 2.]], [[ 3., 3., 3.], [ 3., 3., 3.], [ 3., 3., 3.]]])
So far so good? Here is the average for each of those 3x3 blocks for the first band (just scrunched up the answer)
>>> np.mean(a[0]) # 17.5 overall for band 1 >>> np.mean(a[0,:3,:3]) # 7.0 top left 3x3 block band 0 >>> np.mean(a[1,:3,:3]) # 2.0 top left 3x3 block band 1 >>> np.mean(a[2,:3,:3]) # 3.0 top left 3x3 block band 2
Now you normally work with blocks based on the number of rows and columns rather than 'stride size'. But you can block an array based on the band and the number of rows and columns rather than that arbitrary unexplained 512 block.
here is some code to work with blocks to get you thinks about how to tackle to problem
import numpy as np from numpy.lib.stride_tricks import as_strided def block_array(a, block=(3,3)): """Provide a 2D block view of a 2D array or a slice of a 3D array. No error checking made. Columns and rows outside of the block are truncated. """ r, c = block # 3x3 block default shape = (a.shape[0]/r, a.shape[1]/c) + block strides = (r*a.strides[0], c*a.strides[1]) + a.strides b_a = as_strided(a, shape=shape, strides=strides) return b_a # assign all 3a = np.ones((3,6,6)) # make a 3D array a[0] = np.arange(36).reshape(6,6) # 36 sequential numbers a[1] = 2 # assign all 2 a[2] = 3 block = (3,3) b_0 = block_array(a[0],block) b_1 = block_array(a[1],block) b_2 = block_array(a[2],block)
ps...
I forgot to add, that the simple average through all 3 bands is simply
>>> np.set_printoptions(precision=2,threshold=10,edgeitems=3,linewidth=75,suppress=True) >>> np.mean(a,axis=0,keepdims=False) array([[ 1.67, 2. , 2.33, 2.67, 3. , 3.33], [ 3.67, 4. , 4.33, 4.67, 5. , 5.33], [ 5.67, 6. , 6.33, 6.67, 7. , 7.33], [ 7.67, 8. , 8.33, 8.67, 9. , 9.33], [ 9.67, 10. , 10.33, 10.67, 11. , 11.33], [ 11.67, 12. , 12.33, 12.67, 13. , 13.33]]) >>>
I solved the problem by changing some lines of code... here is the result:
Dan Patterson thanks for the suggestions. I will keep in mind for future elaborations.
Please share your fixed code, so others can benefit. Thanks!