[gdal-dev] ReadAsArray

[Yuta Sato]

On Fri, Apr 3, 2015 at 11:08 PM, Stefan Keller <sfkeller at gmail.com> wrote:

> Hi Yuta and Even
>
> Excuse me when I creep in that discussion - but I'm acutally having
> quite a similar beginners question (I'm rather a vector guy :-) ):
>
> Yuta: Can you post the code fragment here in order to make benchmarks?
>

> I have a GeoTIFF (representing heights from SRTM3) in Mercator CRS of
> the size of a country.
> And I have a position (in lat/lon) and a radius (in m) as input.
> Now, I'd like to read out all grid cells within a radius in order to
> calculate the highest point nearby.
> So I'd like to read out a subset of grid cells/pixels in the most
> efficient way.
>

Hi Stefan,

I think, SRTM3 with country level goes smoothy.
Do you really need to read the data within a radius?


>
> Here's some code snippet:
>   import gdal
>   from gdalconst import *
>   filename = 'my.tif'
>   dataset = gdal.Open(filename, GA_ReadOnly)
>   band = dataset.GetRasterBand(1)
>   scanline = band.ReadRaster( 0, 0, band.XSize, 1,band.XSize, 1,
> band.DataType)
>   ...
>
> Following questions arise:
> 1. Are there alternatives reading out the sub-matrix
>

Alternatively, you could use: http://www.gdal.org/gdallocationinfo.html


> 2. How to calculate the grid cell position in Mercator and lat/lon?


I am reading the data within 3 by 3 pixels square from the central lon/lat
as follows (code adapted from Even):

from osgeo import gdal
from osgeo import osr

ds = gdal.Open(image)

srs = osr.SpatialReference()
srs.ImportFromWkt(ds.GetProjection())

srs_wgs84 = osr.SpatialReference()
srs_wgs84.SetFromUserInput('WGS84')

ct = osr.CoordinateTransformation(srs_wgs84, srs)

for lon, lat in zip(lons, lats):

        (X, Y, _) = ct.TransformPoint(lon, lat)

        geomatrix = ds.GetGeoTransform()
        (success, inv_geometrix) = gdal.InvGeoTransform(geomatrix)
        x = int(inv_geometrix[0] + inv_geometrix[1] * X + inv_geometrix[2]
* Y) - 1
        y = int(inv_geometrix[3] + inv_geometrix[4] * X + inv_geometrix[5]
* Y) - 1

        vals = ds.ReadAsArray(x,y,3,3)

If you get confused with some lines, please let me know. It's using
gdal/osr only though.
I'm sorry if I was mistaken your question.

yuta


> Yours, Stefan
>
>
> 2015-04-03 13:16 GMT+02:00 Yuta Sato <yutaxsato at gmail.com>:
> > Thank you very much "Even Rouault" for making me understood.
> >
> > On Fri, Apr 3, 2015 at 8:13 PM, Even Rouault <even.rouault at spatialys.com
> >
> > wrote:
> >>
> >> Le vendredi 03 avril 2015 12:46:48, Yuta Sato a écrit :
> >> > Dear Even Rouault:
> >> >
> >> > Thank you very much.
> >> > What about setting these parameters "used with .ReadAsArray()",
> though I
> >> > did not know their meanings?
> >> > buf_xsize=None, buf_ysize=None, buf_obj=None
> >>
> >> (Please keep the list CC'ed)
> >>
> >> buf_xsize and buf_ysize are the equivalents of nBufXSize and nBufYSize
> in
> >> GDALRasterBand::RasterIO()
> >>
> http://gdal.org/classGDALRasterBand.html#a75d4af97b3436a4e79d9759eedf89af4
> >> i.e. to do downsampling or upsampling of original data.
> >>
> >> buf_obj can be used to "recycle" an existing numpy array of the
> >> appropriate
> >> size instead of allocating a new one.
> >>
> >> >
> >> >
> >> > On Fri, Apr 3, 2015 at 7:43 PM, Even Rouault
> >> > <even.rouault at spatialys.com>
> >> >
> >> > wrote:
> >> > > Le vendredi 03 avril 2015 12:22:00, Yuta Sato a écrit :
> >> > > > Dear Respected GDAL Developers and Users:
> >> > > >
> >> > > > What parameters should I set beforehand in order to accelerate the
> >> > >
> >> > > reading
> >> > >
> >> > > > of a GeoTiff file?
> >> > > >
> >> > > > I am using as follows:
> >> > > >
> >> > > > data =
> >> > > > src_dataset.GetRasterBand(1).ReadAsArray(xoff,yoff,xsize,ysize)
> >> > > >
> >> > > > Does setting the following parameters accelerate?
> >> > > >
> >> > > > GDAL_CACHEMAX, GDAL_SWATH_SIZE
> >> > > >
> >> > > > I'm using gdal python.
> >> > >
> >> > > Yuta,
> >> > >
> >> > > Increasing GDAL_CACHEMAX might accelerate in case of repeated reads
> on
> >> > > windows
> >> > > that are identical or overlapping already read windows. Or if the
> way
> >> > > you
> >> > > read
> >> > > the raster doesn't follow its block shape : for example if the
> raster
> >> > > is
> >> > > organized by lines/strips and you read by square blocks, or the
> >> > > reverse
> >> > > situation.
> >> > >
> >> > > GDAL_SWATH_SIZE is only used by CreateCopy().
> >> > >
> >> > > Even
> >> > >
> >> > > --
> >> > > Spatialys - Geospatial professional services
> >> > > http://www.spatialys.com
> >>
> >> --
> >> Spatialys - Geospatial professional services
> >> http://www.spatialys.com
> >
> >
> >
> > _______________________________________________
> > gdal-dev mailing list
> > gdal-dev at lists.osgeo.org
> > http://lists.osgeo.org/mailman/listinfo/gdal-dev
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20150404/b5ad6142/attachment-0001.html>