I typically read large raster data using gdal and setting up a moving window of 256*256 pixels (for example).

I recently had to do some processing using arcpy, so I did the same thing, using arcpy.RasterToNumPyArray.

This is incredibly slow compared to using the gdal equivalent. It is almost worth the hassle of exporting the raster to a gdal format and using gdal.

I was wondering why it is so slow. My only thought is that with each call to arcpy.RasterToNumPyArray, you must provide the file path...I imagine arcpy is then opening/closing the file and/or acquiring/releasing locks on each call. Seems like this could waste a certain amount of time.

Is there a more efficient method for iterating through windows of raster data in arcpy?

Here is an example class which will iterate through an ESRI raster in windows:

  class EsriRasterDataset(object):
  Represents an ESRI raster dataset.
  Aims to be somewhat similar to the gdal interface, although provides more
  functions for accessing individual properties of the geotransform
  def __init__(self, path): 

         path: ESRI-style path to the raster (e.g. "C:\\temp\\mygeodatabase.gdb\\myraster")
      self.path = path

  def RasterXSize(self):
      return int(arcpy.GetRasterProperties_management(self.path, "COLUMNCOUNT")[0])

  def RasterYSize(self):
      return int(arcpy.GetRasterProperties_management(self.path, "ROWCOUNT")[0])        

  def GetGeoTransform(self):
      x_res = float(arcpy.GetRasterProperties_management(self.path, "CELLSIZEX")[0])
      y_res = float(arcpy.GetRasterProperties_management(self.path, "CELLSIZEY")[0])

      ulx = float(arcpy.GetRasterProperties_management(self.path, "LEFT")[0])
      uly = float(arcpy.GetRasterProperties_management(self.path, "TOP")[0])

      return [ulx, x_res, 0, uly, 0, y_res]

  def x_res(self):
      x_res = float(arcpy.GetRasterProperties_management(self.path, "CELLSIZEX")[0])
      return x_res

  def y_res(self):
      y_res = float(arcpy.GetRasterProperties_management(self.path, "CELLSIZEY")[0])
      return y_res

  def llx(self):
      return self.ulx()

  def lly(self):
      return float(arcpy.GetRasterProperties_management(self.path, "BOTTOM")[0])

  def ulx(self):
      return float(arcpy.GetRasterProperties_management(self.path, "LEFT")[0])

  def uly(self):
      return float(arcpy.GetRasterProperties_management(self.path, "TOP")[0])

  def GetProjectionRef(self):
      '''TODO: '''

  def ReadAsArray(self, xOff=None, yOff=None, ncols=None, nrows=None):
      '''TODO: '''

  def iter_raster(self, winSize):
      '''generator function to provide an iterator over the dataset'''
      # Bring some properties in scope so we don't have to keep accessing them in a loop
      x_res = self.x_res()
      y_res = self.y_res()
      ysize = self.RasterYSize()
      xsize = self.RasterXSize()

      nxWindows = xsize / winSize
      nyWindows = ysize / winSize

      llx = self.llx()
      lly = self.lly()

      nX = 0        
      # Use a while loop to avoid having to create a potentially huge range() list
      while nX < nxWindows:
        nY = 0
        while nY < nyWindows:
          xOff = (nX * winSize)
          yOff = (nY * winSize)

          xOff_mu = llx + (xOff * x_res)
          yOff_mu = lly + (yOff * y_res)
          origin = arcpy.Point(xOff_mu, yOff_mu)                          
          ncols = winSize
          nrows = winSize
          if winSize > (xsize - xOff):
              ncols = xsize - xOff
          if winSize > (ysize - yOff):
              nrows = ysize - yOff

          yield arcpy.RasterToNumPyArray(self.path, origin, ncols, nrows), ncols, nrows, xOff, yOff
          nX += 1
        nY += 1

The iter_raster function here is more or less the same as I do it when using gdal, although the arcpy interface is more clumsy (it requires the coordinate of the start cell...so presumably has to convert these back to indices. GDAL simply asks for the offset and number of cells...).

Would there be any way that arcpy can simply provide the binary 'scanline' (in the same manner that GDAL will if just calling ReadRaster) rather than converting to a numpy array? This is just a thought. I'm open to anything that improves performance...

(ArcObjects currently seems the only way to go)

  • 3
    It would be helpful to edit your post to include an example script highlighting the functionality you are trying to improve. – Aaron Jun 24 '15 at 16:29
  • 2
    Yes, Esri is slower than GDAL, even though Esri now uses GDAL drivers... there's a whole list of inheritance/objects that the data needs to be converted through before it becomes a numpy array... what rasters are you trying to read? VRT is particularly slow and Mosaic Datasets need to do a lot of processing before you can get any values out... I suspect that it's one of the Esri internal formats or you would just use GDAL. If you want speed then consider ArcObjects, it's much faster (but still not as fast as GDAL). – Michael Stimson Jun 24 '15 at 21:15
  • 1
    @MichaelMiles-Stimson Yes I'm trying to read a simple file geodatabase raster. I want to recreate the raster in PostGIS and so am using GDAL to do the writing, and arcpy to do the reading from the esri raster. I want to provide users a means of doing the conversion in a single step instead of the current method of converting to an intermeadiate format (tif) – James Jun 25 '15 at 7:32
  • 1
    @Aaron added an example – James Jun 25 '15 at 7:40
  • 1
    Another thing to consider is that NumPy arrays are intended to be created once and have many operations run against it, rather than trying to create many arrays that have a single operation done each. You would likely find large performance gains if you can model your workflow around that concept. If you can't, then NumPy isn't the proper tool for your task, as you suspect. – Evil Genius Jun 25 '15 at 11:31