[gdal-dev] flipping raster in python

[Jamie Adams]

Thanks Even,

I had (sort of) come to that conclusion myself, as I started looking into
what blocksize actually meant.  I'll modify to work in row order.

Thanks!

Jamie

On Fri, Mar 26, 2010 at 2:10 PM, Even Rouault
<even.rouault at mines-paris.org>wrote:

> 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
>
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.osgeo.org/pipermail/gdal-dev/attachments/20100326/080d3821/attachment-0001.html