[gdal-dev] GDAL Python Development: Area Selection

[Poovendran Moodley]

Hi Even,

I'm really appreciative of your help! And yes, gt[2] = gt[4] = 0 :) Here's
the code that I've written thus far (for those that are interested --
hopefully I'm not mistaken in my coding!):
def select(filename, x1, y1, x2, y2):
    dataset = gdal.Open( filename, GA_ReadOnly )
    colStart = int(math.ceil((x1 - gt[0])/gt[1]))
    rowStart = int(math.ceil((y1 - gt[3])/gt[5]))
    colEnd = int(math.ceil((x2 - gt[0])/gt[1]))
    rowEnd = int(math.ceil((y2 - gt[3])/gt[5]))
    band = dataset.GetRasterBand(1)
    n = colEnd - colStart + 1
    colStart = colStart - 1 #Because ReadRaster begins at row 0
    for k in range(rowStart - 1, rowEnd):
        scanline = band.ReadRaster(colStart, k, n, 1, n, 1, GDT_Float32) '''
reading 1 row at a time within specified window '''
        tuple_of_floats = struct.unpack('f' * n, scanline)
        print tuple_of_floats

It took me a while to fully understand how ReadRaster works. I read a post
in the FWTools newsgroup:
http://lists.maptools.org/pipermail/fwtools/2005-September/000138.html
The guy there explained the usage quite well. It's a pity that there isn't
more documentation readily available. Thank you again for your assistance
Even! All of the best to you :)

Kind regards
Pooven

On Wed, May 27, 2009 at 9:35 PM, Even Rouault
<even.rouault at mines-paris.org>wrote:

> Pooven,
>
> You're doing things the right way. GDAL doesn't provide any high-level
> function to extract a region from a bounding box expressed in projected
> coordinates. You have to compute the pixel coordinates of the corner of the
> bounding box yourself by using the inverse of the geotransform matrix
> returned by the GetGeoTransform() method of the dataset.
>
> gt = ds.GetGeoTransform()
> x_pixel = (x_projected - gt[0]) / gt[1]
> y_pixel = (y_projected - gt[3]) / gt[5]
>
> (provided that your raster doesn't has any rotation or shearing, that is to
> say that gt[2] = gt[4] = 0)
>
> Best regards,
> Even
>
> Le Wednesday 27 May 2009 12:06:10 Poovendran Moodley, vous avez écrit :
> > Hi there,
> >
> > I want to select only a region of a raster file; that is, suppose I know
> > the coordinates that form a bounding box around a region, I want to
> select
> > just that region and write it to another file.
> >
> > I'm rather new to GDAL and thus far I'm only aware of one way to do this;
> >
> > since I'm using Python:
> > > scanline = band.ReadRaster( 0, 0, band.XSize, 1, band.XSize, 1,
> > > GDT_Float32
> > > <
> http://www.gdal.org/gdal_8h.html#22e22ce0a55036a96f652765793fb7a4f5cbd2f
> > >96abffd9ac061fc0dced5cbba> )
> > >
> > > I believe that I can specify the window of raster data to read using
> the
> >
> > 3rd till 6th input parameters of the above method (with reasonable size
> > reads untill the entire regions I'm interested in is read).
> >
> > Is there some other way this was suppose to be done? The method above
> will
> > require me to convert my coordinates into rows and columns (I think)
> which
> > is fine... but leaves room for error so I'd like to only rely on GDAL
> > routines. Can someone please advise? It could very well be that I'm
> totally
> > mistaken... I hope this isn't such a bad question...
> >
> > Kind regards
> > Pooven
>
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.osgeo.org/pipermail/gdal-dev/attachments/20090529/c8ec325d/attachment.html