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 <[email protected]>: > > > 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
