[gdal-dev] gdallocationinfo.exe multiple access

[Even Rouault]

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