[gdal-dev] GDAL Python Development: Area Selection

Even Rouault even.rouault at mines-paris.org
Fri May 29 13:59:57 EDT 2009


Pooven,

We clearly miss a specific documentation for the GDAL Python bindings 
(contributions welcome !). The OGR ones have docstring help generated from 
C++ doxygen.

Generally, the Python API is close to the C++ one that is documented on the 
website, but there are a few specific methods like ReadRaster() that are 
specific to the binding. There's also the Javadoc of the Java bindings 
(http://gdal.org/java/) that can be somehow useful, as they share a lot of 
things in common, but expect a few differences too.

A few other potential helpfull pointers :
* Someone brought my attention to the following site, 
http://www.gis.usu.edu/~chrisg/python/, that has very interesting slides and 
exercices about how to use GDAL/OGR in Python.
* The various utilities written in Python, can also serve as exemple : 
http://trac.osgeo.org/gdal/browser/trunk/gdal/swig/python/scripts
http://trac.osgeo.org/gdal/browser/trunk/gdal/swig/python/samples

(Ah, and Frank Warmerdam, "the guy there" answering in the post you refer to, 
is the main author of GDAL ;-) )

Even

Le Friday 29 May 2009 08:43:08 Poovendran Moodley, vous avez écrit :
> 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




More information about the gdal-dev mailing list