[gdal-dev] ReadAsArray

[Stefan Keller]

Hi Yuta

Thanks! I'll try this.
I'm still wondering which approach is faster: With gdal or rasterio...

Cheers, Stefan


2015-04-03 17:46 GMT+02:00 Yuta Sato <yutaxsato at gmail.com>:
>
>
> 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
>
>