[gdal-dev] GDAL Python Development: Area Selection

[Poovendran Moodley]

Hi Even,

Thank you for the additional resources! hehe and thanks to Frank for making
GDAL available :)

Kind regards
Pooven

On Fri, May 29, 2009 at 7:59 PM, Even Rouault
<even.rouault at mines-paris.org>wrote:

> 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/<http://www.gis.usu.edu/%7Echrisg/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
>
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.osgeo.org/pipermail/gdal-dev/attachments/20090602/92ea08a6/attachment-0001.html