[gdal-dev] flipping raster in python

[Even Rouault]

Jamie,

I'm not sure why it slowdowns after the first 2 bands, but I can say your way 
of proceeding is very inefficient and you can achieve the same results by 
processing by rows instead of columns. As your input file is scanline 
oriented (block size is 20000x1 - which is pretty standard), reading a column 
will cause the whole file to be read at each time !

Try instead the following version which read line by line and horizontally 
flip the buffer : it should complete in a few minutes.

for band in range(1, indataset.RasterCount + 1):
  ## Print the current band
  print 'Band %s' % band

  inband = indataset.GetRasterBand(band)
  outband = outdataset.GetRasterBand(band)

  count = 0
  print count

  for i in range(inband.YSize):

    ## print the current interval of 1k lines
    new_count = int(i / 1000)
    if new_count > count:
      print new_count
      count = new_count

    inline = inband.ReadAsArray(0, i, inband.XSize, 1)
    flipline = numpy.fliplr(inline)
    outband.WriteArray(flipline, 0, i)

    inline = None

  inband = None
  outband = None

Best regards,

Even

Le Friday 26 March 2010 21:22:54 Jamie Adams, vous avez écrit :
> I've written some simple Python code to flip a raster in relation to the
> y-axis.  The raster is 20000x19459 and has 4 bands of type Byte, and is
> written east to west (I have no idea why).  The script proceeds normally
> for the first 2 bands, but slows way down after starting band 3.  I let it
> run overnight, and it had processed less than 1000 columnsin band 3,
> whereas the previous pace for bands 1&2 was around 1000 columns every 10
> seconds.  Here is the relevant code with a few debugging statements
> followed by the gdalinfo output.  The only thing I'm noticing is that the
> input image is INTERLEAVED=BAND whereas the output is INTERLEAVED=PIXEL. 
> Any help is appreciated.
>
> ---------------------------------------------------------------------------
>--- for band in range(1, indataset.RasterCount + 1):
>   ## Print the current band
>   print 'Band %s' % band
>
>   inband = indataset.GetRasterBand(band)
>   outband = outdataset.GetRasterBand(band)
>
>   count = 0
>   print count
>
>   for i in range(inband.XSize):
>
>     ## print the current interval of 1k lines
>     new_count = int(i / 1000)
>     if new_count > count:
>       print new_count
>       count = new_count
>
>     inline = inband.ReadAsArray(i, 0, 1, inband.YSize)
>
>     outband.WriteArray(inline, inband.XSize - i - 1, 0)
>
>     inline = None
>
>   inband = None
>   outband = None
>
> ---------------------------------------------------------------------------
>--- Driver: GTiff/GeoTIFF
> Files: test.tif
> Size is 20000, 19459
> Coordinate System is:
> GEOGCS["WGS 84",
>     DATUM["WGS_1984",
>         SPHEROID["WGS 84",6378137,298.257223563,
>             AUTHORITY["EPSG","7030"]],
>         AUTHORITY["EPSG","6326"]],
>     PRIMEM["Greenwich",0],
>     UNIT["degree",0.0174532925199433],
>     AUTHORITY["EPSG","4326"]]
> Origin = (-175.536864949200321,33.666235092496755)
> Pixel Size = (-0.000062579160220,-0.000054370391310)
> Metadata:
>   AREA_OR_POINT=Area
> Image Structure Metadata:
>   INTERLEAVE=BAND
> Corner Coordinates:
> Upper Left  (-175.5368649,  33.6662351) (175d32'12.71"W, 33d39'58.45"N)
> Lower Left  (-175.5368649,  32.6082416) (175d32'12.71"W, 32d36'29.67"N)
> Upper Right (-176.7884482,  33.6662351) (176d47'18.41"W, 33d39'58.45"N)
> Lower Right (-176.7884482,  32.6082416) (176d47'18.41"W, 32d36'29.67"N)
> Center      (-176.1626566,  33.1372384) (176d 9'45.56"W, 33d 8'14.06"N)
> Band 1 Block=20000x1 Type=Byte, ColorInterp=Red
>   Mask Flags: PER_DATASET ALPHA
> Band 2 Block=20000x1 Type=Byte, ColorInterp=Green
>   Mask Flags: PER_DATASET ALPHA
> Band 3 Block=20000x1 Type=Byte, ColorInterp=Blue
>   Mask Flags: PER_DATASET ALPHA
> Band 4 Block=20000x1 Type=Byte, ColorInterp=Alpha