title: gdal ReadAsarray for vrt extremely slow
body: <p>I am trying to simply subtract two rasters and save the result in another raster. One of the input images is a tif-file, the other is a vrt-file. (Output is tif)</p>
  
  <p>The files are very big, so I open them, divide them into tiles and run through each of them and then subtracting. <strong>The problem is that it is extremely slow!</strong></p>
  
  <pre><code>import gdal
  import numpy as np
  
  rA = gdal.Open(raFileName)
  rB = gdal.Open(rbFileName)
  
  nodataA = rA.GetRasterBand(1).GetNoDataValue()
  nodataB = rB.GetRasterBand(1).GetNoDataValue()
  
  
  raster = gdal.GetDriverByName('GTiff').Create(outputFileName, ncols, nrows, 1 ,gdal.GDT_Float32,['COMPRESS=LZW','BigTIFF=YES'])
  
  # tile size
  trows = 5000
  tcols = 5000
  # number of tiles in output file   (ceil)
  ntrows = (nrows-1)/trows+1
  ntcols = (ncols-1)/tcols+1
  
  # tiling because of memory problems
  for r in range(ntrows):
      for c in range(ntcols):
          # number of rows/cols for tile r (in case of edge)
          rtrows = min(trows,nrows-r*trows)
          ctcols = min(tcols,ncols-c*tcols)
  
          # the data from the input files
          rA_arr = rA.GetRasterBand(1).ReadAsArray(c*tcols,r*trows,ctcols,rtrows)
          rB_arr = rB.GetRasterBand(1).ReadAsArray(c*tcols,r*trows,ctcols,rtrows)
  
          mask = np.logical_or(rA_arr==nodataA,rB_arr==nodataB)
  
          subtraction = rA_arr-rB_arr
          subtraction[mask] = nodata
  
          # writing this tile to the output raster
          rasterBand.WriteArray(subtraction,c*tcols,r*trows)            
          raster.FlushCache()
  </code></pre>
  
  <p>The rasters I am currently trying to subtract have 16*16 tiles (of 5000*5000 pxls) and after about two hours, it has only gone through 3 rows!</p>
  
  <p><strong>Is there any way to increase performance??</strong></p>