Hi Even,<br><br>Thank you for the additional resources! hehe and thanks to Frank for making GDAL available :)<br><br clear="all">Kind regards<br>Pooven<br><br><div class="gmail_quote">On Fri, May 29, 2009 at 7:59 PM, Even Rouault <span dir="ltr"><<a href="mailto:even.rouault@mines-paris.org">even.rouault@mines-paris.org</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">Pooven,<br>
<br>
We clearly miss a specific documentation for the GDAL Python bindings<br>
(contributions welcome !). The OGR ones have docstring help generated from<br>
C++ doxygen.<br>
<br>
Generally, the Python API is close to the C++ one that is documented on the<br>
website, but there are a few specific methods like ReadRaster() that are<br>
specific to the binding. There's also the Javadoc of the Java bindings<br>
(<a href="http://gdal.org/java/" target="_blank">http://gdal.org/java/</a>) that can be somehow useful, as they share a lot of<br>
things in common, but expect a few differences too.<br>
<br>
A few other potential helpfull pointers :<br>
* Someone brought my attention to the following site,<br>
<a href="http://www.gis.usu.edu/%7Echrisg/python/" target="_blank">http://www.gis.usu.edu/~chrisg/python/</a>, that has very interesting slides and<br>
exercices about how to use GDAL/OGR in Python.<br>
* The various utilities written in Python, can also serve as exemple :<br>
<a href="http://trac.osgeo.org/gdal/browser/trunk/gdal/swig/python/scripts" target="_blank">http://trac.osgeo.org/gdal/browser/trunk/gdal/swig/python/scripts</a><br>
<a href="http://trac.osgeo.org/gdal/browser/trunk/gdal/swig/python/samples" target="_blank">http://trac.osgeo.org/gdal/browser/trunk/gdal/swig/python/samples</a><br>
<br>
(Ah, and Frank Warmerdam, "the guy there" answering in the post you refer to,<br>
is the main author of GDAL ;-) )<br>
<br>
Even<br>
<br>
Le Friday 29 May 2009 08:43:08 Poovendran Moodley, vous avez écrit :<br>
<div><div></div><div class="h5">> Hi Even,<br>
><br>
> I'm really appreciative of your help! And yes, gt[2] = gt[4] = 0 :) Here's<br>
> the code that I've written thus far (for those that are interested --<br>
> hopefully I'm not mistaken in my coding!):<br>
> def select(filename, x1, y1, x2, y2):<br>
> dataset = gdal.Open( filename, GA_ReadOnly )<br>
> colStart = int(math.ceil((x1 - gt[0])/gt[1]))<br>
> rowStart = int(math.ceil((y1 - gt[3])/gt[5]))<br>
> colEnd = int(math.ceil((x2 - gt[0])/gt[1]))<br>
> rowEnd = int(math.ceil((y2 - gt[3])/gt[5]))<br>
> band = dataset.GetRasterBand(1)<br>
> n = colEnd - colStart + 1<br>
> colStart = colStart - 1 #Because ReadRaster begins at row 0<br>
> for k in range(rowStart - 1, rowEnd):<br>
> scanline = band.ReadRaster(colStart, k, n, 1, n, 1, GDT_Float32)<br>
> ''' reading 1 row at a time within specified window '''<br>
> tuple_of_floats = struct.unpack('f' * n, scanline)<br>
> print tuple_of_floats<br>
><br>
> It took me a while to fully understand how ReadRaster works. I read a post<br>
> in the FWTools newsgroup:<br>
> <a href="http://lists.maptools.org/pipermail/fwtools/2005-September/000138.html" target="_blank">http://lists.maptools.org/pipermail/fwtools/2005-September/000138.html</a><br>
> The guy there explained the usage quite well. It's a pity that there isn't<br>
> more documentation readily available. Thank you again for your assistance<br>
> Even! All of the best to you :)<br>
><br>
> Kind regards<br>
> Pooven<br>
><br>
> On Wed, May 27, 2009 at 9:35 PM, Even Rouault<br>
><br>
> <<a href="mailto:even.rouault@mines-paris.org">even.rouault@mines-paris.org</a>>wrote:<br>
> > Pooven,<br>
> ><br>
> > You're doing things the right way. GDAL doesn't provide any high-level<br>
> > function to extract a region from a bounding box expressed in projected<br>
> > coordinates. You have to compute the pixel coordinates of the corner of<br>
> > the bounding box yourself by using the inverse of the geotransform matrix<br>
> > returned by the GetGeoTransform() method of the dataset.<br>
> ><br>
> > gt = ds.GetGeoTransform()<br>
> > x_pixel = (x_projected - gt[0]) / gt[1]<br>
> > y_pixel = (y_projected - gt[3]) / gt[5]<br>
> ><br>
> > (provided that your raster doesn't has any rotation or shearing, that is<br>
> > to say that gt[2] = gt[4] = 0)<br>
> ><br>
> > Best regards,<br>
> > Even<br>
> ><br>
> > Le Wednesday 27 May 2009 12:06:10 Poovendran Moodley, vous avez écrit :<br>
> > > Hi there,<br>
> > ><br>
> > > I want to select only a region of a raster file; that is, suppose I<br>
> > > know the coordinates that form a bounding box around a region, I want<br>
> > > to<br>
> ><br>
> > select<br>
> ><br>
> > > just that region and write it to another file.<br>
> > ><br>
> > > I'm rather new to GDAL and thus far I'm only aware of one way to do<br>
> > > this;<br>
> > ><br>
> > > since I'm using Python:<br>
> > > > scanline = band.ReadRaster( 0, 0, band.XSize, 1, band.XSize, 1,<br>
> > > > GDT_Float32<br>
> > > > <<br>
> ><br>
> > <a href="http://www.gdal.org/gdal_8h.html#22e22ce0a55036a96f652765793fb7a4f5cbd2f" target="_blank">http://www.gdal.org/gdal_8h.html#22e22ce0a55036a96f652765793fb7a4f5cbd2f</a><br>
> ><br>
> > > >96abffd9ac061fc0dced5cbba> )<br>
> > > ><br>
> > > > I believe that I can specify the window of raster data to read using<br>
> ><br>
> > the<br>
> ><br>
> > > 3rd till 6th input parameters of the above method (with reasonable size<br>
> > > reads untill the entire regions I'm interested in is read).<br>
> > ><br>
> > > Is there some other way this was suppose to be done? The method above<br>
> ><br>
> > will<br>
> ><br>
> > > require me to convert my coordinates into rows and columns (I think)<br>
> ><br>
> > which<br>
> ><br>
> > > is fine... but leaves room for error so I'd like to only rely on GDAL<br>
> > > routines. Can someone please advise? It could very well be that I'm<br>
> ><br>
> > totally<br>
> ><br>
> > > mistaken... I hope this isn't such a bad question...<br>
> > ><br>
> > > Kind regards<br>
> > > Pooven<br>
<br>
<br>
</div></div></blockquote></div><br>