On Fri, Apr 3, 2015 at 11:08 PM, Stefan Keller <[email protected]> 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 <[email protected]>: > > Thank you very much "Even Rouault" for making me understood. > > > > On Fri, Apr 3, 2015 at 8:13 PM, Even Rouault <[email protected] > > > > 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 > >> > <[email protected]> > >> > > >> > 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 > > [email protected] > > http://lists.osgeo.org/mailman/listinfo/gdal-dev >
_______________________________________________ gdal-dev mailing list [email protected] http://lists.osgeo.org/mailman/listinfo/gdal-dev
