<div dir="ltr"><br><div class="gmail_extra"><br><div class="gmail_quote">On Fri, Apr 3, 2015 at 11:08 PM, Stefan Keller <span dir="ltr"><<a href="mailto:sfkeller@gmail.com" target="_blank">sfkeller@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">Hi Yuta and Even<br>
<br>
Excuse me when I creep in that discussion - but I'm acutally having<br>
quite a similar beginners question (I'm rather a vector guy :-) ):<br>
<br>
Yuta: Can you post the code fragment here in order to make benchmarks?<br></blockquote><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">
<br>
I have a GeoTIFF (representing heights from SRTM3) in Mercator CRS of<br>
the size of a country.<br>
And I have a position (in lat/lon) and a radius (in m) as input.<br>
Now, I'd like to read out all grid cells within a radius in order to<br>
calculate the highest point nearby.<br>
So I'd like to read out a subset of grid cells/pixels in the most efficient way.<br></blockquote><div><br>Hi Stefan,<br><br>I think, SRTM3 with country level goes smoothy. <br>Do you really need to read the data within a radius?<br> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">
<br>
Here's some code snippet:<br>
  import gdal<br>
  from gdalconst import *<br>
  filename = 'my.tif'<br>
  dataset = gdal.Open(filename, GA_ReadOnly)<br>
  band = dataset.GetRasterBand(1)<br>
  scanline = band.ReadRaster( 0, 0, band.XSize, 1,band.XSize, 1, band.DataType)<br>
  ...<br>
<br>
Following questions arise:<br>
1. Are there alternatives reading out the sub-matrix<br></blockquote><div><br>Alternatively, you could use: <a href="http://www.gdal.org/gdallocationinfo.html">http://www.gdal.org/gdallocationinfo.html</a><br> <br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">
2. How to calculate the grid cell position in Mercator and lat/lon?</blockquote><div><br>I
 am reading the data within 3 by 3 pixels square from the central 
lon/lat as follows (code adapted from Even):<br><br>from osgeo import gdal<br>
from osgeo import osr<br>
<br>
ds = gdal.Open(image)<br>
<br>
srs = osr.SpatialReference()<br>
srs.ImportFromWkt(ds.GetProjection())<br>
<br>
srs_wgs84 = osr.SpatialReference()<br>
srs_wgs84.SetFromUserInput('WGS84')<br>
<br>
ct = osr.CoordinateTransformation(srs_wgs84, srs)   <br>
<br>
for lon, lat in zip(lons, lats):<br>
<br>
        (X, Y, _) = ct.TransformPoint(lon, lat)<br>
<br>
        geomatrix = ds.GetGeoTransform()<br>
        (success, inv_geometrix) = gdal.InvGeoTransform(geomatrix)<br>
        x = int(inv_geometrix[0] + inv_geometrix[1] * X + inv_geometrix[2] * Y) - 1<br>
        y = int(inv_geometrix[3] + inv_geometrix[4] * X + inv_geometrix[5] * Y) - 1<br>
<br>
        vals = ds.ReadAsArray(x,y,3,3)<br><div><br>If you get confused with some lines, please let me know. It's using gdal/osr only though.<br></div><div>I'm sorry if I was mistaken your question.<br><br></div><div>yuta <br><br></div></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">
<br>
Yours, Stefan<br>
<div><div class="h5"><br>
<br>
2015-04-03 13:16 GMT+02:00 Yuta Sato <<a href="mailto:yutaxsato@gmail.com">yutaxsato@gmail.com</a>>:<br>
> Thank you very much "Even Rouault" for making me understood.<br>
><br>
> On Fri, Apr 3, 2015 at 8:13 PM, Even Rouault <<a href="mailto:even.rouault@spatialys.com">even.rouault@spatialys.com</a>><br>
> wrote:<br>
>><br>
>> Le vendredi 03 avril 2015 12:46:48, Yuta Sato a écrit :<br>
>> > Dear Even Rouault:<br>
>> ><br>
>> > Thank you very much.<br>
>> > What about setting these parameters "used with .ReadAsArray()", though I<br>
>> > did not know their meanings?<br>
>> > buf_xsize=None, buf_ysize=None, buf_obj=None<br>
>><br>
>> (Please keep the list CC'ed)<br>
>><br>
>> buf_xsize and buf_ysize are the equivalents of nBufXSize and nBufYSize in<br>
>> GDALRasterBand::RasterIO()<br>
>> <a href="http://gdal.org/classGDALRasterBand.html#a75d4af97b3436a4e79d9759eedf89af4" target="_blank">http://gdal.org/classGDALRasterBand.html#a75d4af97b3436a4e79d9759eedf89af4</a><br>
>> i.e. to do downsampling or upsampling of original data.<br>
>><br>
>> buf_obj can be used to "recycle" an existing numpy array of the<br>
>> appropriate<br>
>> size instead of allocating a new one.<br>
>><br>
>> ><br>
>> ><br>
>> > On Fri, Apr 3, 2015 at 7:43 PM, Even Rouault<br>
>> > <<a href="mailto:even.rouault@spatialys.com">even.rouault@spatialys.com</a>><br>
>> ><br>
>> > wrote:<br>
>> > > Le vendredi 03 avril 2015 12:22:00, Yuta Sato a écrit :<br>
>> > > > Dear Respected GDAL Developers and Users:<br>
>> > > ><br>
>> > > > What parameters should I set beforehand in order to accelerate the<br>
>> > ><br>
>> > > reading<br>
>> > ><br>
>> > > > of a GeoTiff file?<br>
>> > > ><br>
>> > > > I am using as follows:<br>
>> > > ><br>
>> > > > data =<br>
>> > > > src_dataset.GetRasterBand(1).ReadAsArray(xoff,yoff,xsize,ysize)<br>
>> > > ><br>
>> > > > Does setting the following parameters accelerate?<br>
>> > > ><br>
>> > > > GDAL_CACHEMAX, GDAL_SWATH_SIZE<br>
>> > > ><br>
>> > > > I'm using gdal python.<br>
>> > ><br>
>> > > Yuta,<br>
>> > ><br>
>> > > Increasing GDAL_CACHEMAX might accelerate in case of repeated reads on<br>
>> > > windows<br>
>> > > that are identical or overlapping already read windows. Or if the way<br>
>> > > you<br>
>> > > read<br>
>> > > the raster doesn't follow its block shape : for example if the raster<br>
>> > > is<br>
>> > > organized by lines/strips and you read by square blocks, or the<br>
>> > > reverse<br>
>> > > situation.<br>
>> > ><br>
>> > > GDAL_SWATH_SIZE is only used by CreateCopy().<br>
>> > ><br>
>> > > Even<br>
>> > ><br>
>> > > --<br>
>> > > Spatialys - Geospatial professional services<br>
>> > > <a href="http://www.spatialys.com" target="_blank">http://www.spatialys.com</a><br>
>><br>
>> --<br>
>> Spatialys - Geospatial professional services<br>
>> <a href="http://www.spatialys.com" target="_blank">http://www.spatialys.com</a><br>
><br>
><br>
><br>
</div></div>> _______________________________________________<br>
> gdal-dev mailing list<br>
> <a href="mailto:gdal-dev@lists.osgeo.org">gdal-dev@lists.osgeo.org</a><br>
> <a href="http://lists.osgeo.org/mailman/listinfo/gdal-dev" target="_blank">http://lists.osgeo.org/mailman/listinfo/gdal-dev</a><br>
</blockquote></div><br></div></div>