[gdal-dev] gdallocationinfo.exe multiple access

[Even Rouault]

Le samedi 28 février 2015 19:03:21, Yuta Sato a écrit :
> Dear Even Rouault:
> Thanks for full Python script. I understood this.
> 
> However, since my lons and lats are thousands, it is taking days to extract
> the values. So, my question is how can I sped up it. Is it possible to call
> several tuples of lons and lats simultaneously? For loop is doing one by
> one right?

Days for just thousands of values ? I hope not...
If you suffer from bad performance, you can try to sort the values by 
decreasing latitudes (and for equal latitude, by increasing longitude), so as 
to be nice with the typical organization of raster files and with the GDAL 
block cache. Alternatively/complementarily you can increase the maximum size 
of the GDAL block cache. The default value is 40 MB. If you call 
gdal.SetConfigOption('GDAL_CACHEMAX', '500'), you'll increase it to 500 MB.
Otherwise if the dataset can entirely fit into memory, you can read it fully 
with ds.ReadAsArray(), and read directly from the returned array with the 
(x,y) computed as below.
Parallel python could be a solution if you really need parallized processing. 
But I guess the complication would only be worth if you must deal with 
millions or billions of coordinate pairs.

> 
> On Sun, Mar 1, 2015 at 2:42 AM, Even Rouault <even.rouault at spatialys.com>
> 
> wrote:
> > Le samedi 28 février 2015 18:29:01, Yuta Sato a écrit :
> > > Dear Even Rouault:
> > > 
> > > Thanks for you quick response.
> > > Unfortunately, I do not know how to pipe using python script.
> > > I just know to use as follows:
> > > 
> > > for x, y in zip(lons, lats):
> > >         vals = subprocess.check_output([gdallocationinfo, '-valonly',
> > > 
> > > '-l_srs', 'wgs84', image, str(x), str(y)])
> > > 
> > > The input lons and lats are in the text file.
> > 
> > As I suggested, you'd better extract the gist of
> > http://svn.osgeo.org/gdal/trunk/gdal/swig/python/samples/val_at_coord.py
> > 
> > So something like
> > 
> > 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)
> > 
> >         y = int(inv_geometrix[3] + inv_geometrix[4] * X +
> >         inv_geometrix[5]
> > 
> > * Y)
> > 
> >         vals = ds.ReadAsArray(x,y,1,1)
> > > 
> > > Can you help me?
> > > <https://plus.google.com/u/0/106381021255842353792?prsrc=4>
> > > 
> > > On Sun, Mar 1, 2015 at 2:05 AM, Even Rouault
> > > <even.rouault at spatialys.com
> > > 
> > > wrote:
> > > > Le samedi 28 février 2015 17:51:47, Yuta Sato a écrit :
> > > > > I have a list of longitudes and latitudes for which I extract the
> > 
> > pixel
> > 
> > > > > values using gdallocationinfo.exe one by one for each pairs of
> > > > > longitudes and latitudes using for loop.
> > > > > 
> > > > > How I simultaneously get pixel values for all the list of
> > > > > longitudes and latitudes?
> > > > > My list is large, so I want faster retrieval.
> > > > > I am using python and outputting the pixel values into text file.
> > > > > Is it possible?
> > > > 
> > > > Yuka,
> > > > 
> > > > you can pipe into gdallocationinfo several tuples of coordinates:
> > > > 
> > > > $ printf "0 0\n0 1\n" | gdallocationinfo byte.tif
> > > > 
> > > > Report:
> > > >   Location: (0P,0L)
> > > >   
> > > >   Band 1:
> > > >     Value: 107
> > > > 
> > > > Report:
> > > >   Location: (0P,1L)
> > > >   
> > > >   Band 1:
> > > >     Value: 115
> > > > 
> > > > But I think it would be easier/faster to just do that fully in
> > > > Python.
> > 
> > > > The following script has likely everything you need to do that:
> > http://svn.osgeo.org/gdal/trunk/gdal/swig/python/samples/val_at_coord.py
> > 
> > > > Even
> > > > 
> > > > 
> > > > --
> > > > Spatialys - Geospatial professional services
> > > > http://www.spatialys.com
> > 
> > --
> > Spatialys - Geospatial professional services
> > http://www.spatialys.com

-- 
Spatialys - Geospatial professional services
http://www.spatialys.com