Hi Even,<br><br>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!):<br>
<span style="font-family: courier new,monospace;">def select(filename, x1, y1, x2, y2):</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> dataset = gdal.Open( filename, GA_ReadOnly )</span><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;"> colStart = int(math.ceil((x1 - gt[0])/gt[1]))</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> rowStart = int(math.ceil((y1 - gt[3])/gt[5]))</span><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;"> colEnd = int(math.ceil((x2 - gt[0])/gt[1]))</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> rowEnd = int(math.ceil((y2 - gt[3])/gt[5]))</span><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;"> band = dataset.GetRasterBand(1)</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> n = colEnd - colStart + 1</span><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;"> colStart = colStart - 1</span> #Because ReadRaster begins at row 0<br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> for k in range(rowStart - 1, rowEnd):</span><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;"> scanline = band.ReadRaster(colStart, k, n, 1, n, 1, GDT_Float32) ''' reading 1 row at a time within specified window '''</span><br style="font-family: courier new,monospace;">
<span style="font-family: courier new,monospace;"> tuple_of_floats = struct.unpack('f' * n, scanline)</span><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;"> print tuple_of_floats</span><br>
<br clear="all">It took me a while to fully understand how ReadRaster works. I read a post in the FWTools newsgroup: <a href="http://lists.maptools.org/pipermail/fwtools/2005-September/000138.html">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 more documentation readily available. Thank you again for your assistance Even! All of the best to you :)<br><br>Kind regards<br>Pooven<br>
<br><div class="gmail_quote">On Wed, May 27, 2009 at 9:35 PM, Even Rouault <span dir="ltr"><<a href="mailto:even.rouault@mines-paris.org" target="_blank">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>
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 the<br>
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 to<br>
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>
<div>> Hi there,<br>
><br>
> I want to select only a region of a raster file; that is, suppose I know<br>
> the coordinates that form a bounding box around a region, I want to select<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 this;<br>
><br>
> since I'm using Python:<br>
> > scanline = band.ReadRaster( 0, 0, band.XSize, 1, band.XSize, 1,<br>
> > GDT_Float32<br>
</div>> > <<a href="http://www.gdal.org/gdal_8h.html#22e22ce0a55036a96f652765793fb7a4f5cbd2f" target="_blank">http://www.gdal.org/gdal_8h.html#22e22ce0a55036a96f652765793fb7a4f5cbd2f</a><br>
> >96abffd9ac061fc0dced5cbba> )<br>
<div><div></div><div>> ><br>
> > I believe that I can specify the window of raster data to read using 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 will<br>
> require me to convert my coordinates into rows and columns (I think) which<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 totally<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>