<div dir="ltr"><div>Yuta, Even:</div><div><br></div><div>Rasterio's rio-sample command</div><div><br></div>  <a href="https://github.com/mapbox/rasterio/blob/master/docs/cli.rst#sample">https://github.com/mapbox/rasterio/blob/master/docs/cli.rst#sample</a> <br><div>  <a href="https://github.com/mapbox/rasterio/blob/master/rasterio/rio/sample.py">https://github.com/mapbox/rasterio/blob/master/rasterio/rio/sample.py</a><br></div><div><br></div><div>can do any number of x, y pairs and only loads image data once. I think it's a pretty good example of how to do this. And it's all GDAL under the hood (but not GDAL's own Python bindings), so you get the same answer as you would from GDAL.</div><div><br></div></div><div class="gmail_extra"><br><div class="gmail_quote">On Sat, Feb 28, 2015 at 11:11 AM, Even Rouault <span dir="ltr"><<a href="mailto:even.rouault@spatialys.com" target="_blank">even.rouault@spatialys.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">Le samedi 28 février 2015 19:03:21, Yuta Sato a écrit :<br>
> Dear Even Rouault:<br>
> Thanks for full Python script. I understood this.<br>
><br>
> However, since my lons and lats are thousands, it is taking days to extract<br>
> the values. So, my question is how can I sped up it. Is it possible to call<br>
> several tuples of lons and lats simultaneously? For loop is doing one by<br>
> one right?<br>
<br>
Days for just thousands of values ? I hope not...<br>
If you suffer from bad performance, you can try to sort the values by<br>
decreasing latitudes (and for equal latitude, by increasing longitude), so as<br>
to be nice with the typical organization of raster files and with the GDAL<br>
block cache. Alternatively/complementarily you can increase the maximum size<br>
of the GDAL block cache. The default value is 40 MB. If you call<br>
gdal.SetConfigOption('GDAL_CACHEMAX', '500'), you'll increase it to 500 MB.<br>
Otherwise if the dataset can entirely fit into memory, you can read it fully<br>
with ds.ReadAsArray(), and read directly from the returned array with the<br>
(x,y) computed as below.<br>
Parallel python could be a solution if you really need parallized processing.<br>
But I guess the complication would only be worth if you must deal with<br>
millions or billions of coordinate pairs.<br>
<br>
><br>
> On Sun, Mar 1, 2015 at 2:42 AM, Even Rouault <<a href="mailto:even.rouault@spatialys.com">even.rouault@spatialys.com</a>><br>
<div class="HOEnZb"><div class="h5">><br>
> wrote:<br>
> > Le samedi 28 février 2015 18:29:01, Yuta Sato a écrit :<br>
> > > Dear Even Rouault:<br>
> > ><br>
> > > Thanks for you quick response.<br>
> > > Unfortunately, I do not know how to pipe using python script.<br>
> > > I just know to use as follows:<br>
> > ><br>
> > > for x, y in zip(lons, lats):<br>
> > >         vals = subprocess.check_output([gdallocationinfo, '-valonly',<br>
> > ><br>
> > > '-l_srs', 'wgs84', image, str(x), str(y)])<br>
> > ><br>
> > > The input lons and lats are in the text file.<br>
> ><br>
> > As I suggested, you'd better extract the gist of<br>
> > <a href="http://svn.osgeo.org/gdal/trunk/gdal/swig/python/samples/val_at_coord.py" target="_blank">http://svn.osgeo.org/gdal/trunk/gdal/swig/python/samples/val_at_coord.py</a><br>
> ><br>
> > So something like<br>
> ><br>
> > from osgeo import gdal<br>
> > from osgeo import osr<br>
> ><br>
> > ds = gdal.Open( image )<br>
> ><br>
> > srs = osr.SpatialReference()<br>
> > srs.ImportFromWkt(ds.GetProjection())<br>
> ><br>
> > srs_wgs84 = osr.SpatialReference()<br>
> > srs_wgs84.SetFromUserInput('WGS84')<br>
> ><br>
> > ct = osr.CoordinateTransformation(srs_wgs84, srs)<br>
> ><br>
> > for lon, lat in zip(lons, lats):<br>
> >         (X, Y, _) = ct.TransformPoint(lon, lat)<br>
> ><br>
> >         geomatrix = ds.GetGeoTransform()<br>
> >         (success, inv_geometrix) = gdal.InvGeoTransform(geomatrix)<br>
> >         x = int(inv_geometrix[0] + inv_geometrix[1] * X +<br>
> >         inv_geometrix[2]<br>
> ><br>
> > * Y)<br>
> ><br>
> >         y = int(inv_geometrix[3] + inv_geometrix[4] * X +<br>
> >         inv_geometrix[5]<br>
> ><br>
> > * Y)<br>
> ><br>
> >         vals = ds.ReadAsArray(x,y,1,1)<br>
> > ><br>
> > > Can you help me?<br>
> > > <<a href="https://plus.google.com/u/0/106381021255842353792?prsrc=4" target="_blank">https://plus.google.com/u/0/106381021255842353792?prsrc=4</a>><br>
> > ><br>
> > > On Sun, Mar 1, 2015 at 2:05 AM, Even Rouault<br>
> > > <<a href="mailto:even.rouault@spatialys.com">even.rouault@spatialys.com</a><br>
> > ><br>
> > > wrote:<br>
> > > > Le samedi 28 février 2015 17:51:47, Yuta Sato a écrit :<br>
> > > > > I have a list of longitudes and latitudes for which I extract the<br>
> ><br>
> > pixel<br>
> ><br>
> > > > > values using gdallocationinfo.exe one by one for each pairs of<br>
> > > > > longitudes and latitudes using for loop.<br>
> > > > ><br>
> > > > > How I simultaneously get pixel values for all the list of<br>
> > > > > longitudes and latitudes?<br>
> > > > > My list is large, so I want faster retrieval.<br>
> > > > > I am using python and outputting the pixel values into text file.<br>
> > > > > Is it possible?<br>
> > > ><br>
> > > > Yuka,<br>
> > > ><br>
> > > > you can pipe into gdallocationinfo several tuples of coordinates:<br>
> > > ><br>
> > > > $ printf "0 0\n0 1\n" | gdallocationinfo byte.tif<br>
> > > ><br>
> > > > Report:<br>
> > > >   Location: (0P,0L)<br>
> > > ><br>
> > > >   Band 1:<br>
> > > >     Value: 107<br>
> > > ><br>
> > > > Report:<br>
> > > >   Location: (0P,1L)<br>
> > > ><br>
> > > >   Band 1:<br>
> > > >     Value: 115<br>
> > > ><br>
> > > > But I think it would be easier/faster to just do that fully in<br>
> > > > Python.<br>
> ><br>
> > > > The following script has likely everything you need to do that:<br>
> > <a href="http://svn.osgeo.org/gdal/trunk/gdal/swig/python/samples/val_at_coord.py" target="_blank">http://svn.osgeo.org/gdal/trunk/gdal/swig/python/samples/val_at_coord.py</a><br>
> ><br>
> > > > Even<br>
> > > ><br>
> > > ><br>
> > > > --<br>
> > > > Spatialys - Geospatial professional services<br>
> > > > <a href="http://www.spatialys.com" target="_blank">http://www.spatialys.com</a><br>
> ><br>
> > --<br>
> > Spatialys - Geospatial professional services<br>
> > <a href="http://www.spatialys.com" target="_blank">http://www.spatialys.com</a><br>
<br>
--<br>
Spatialys - Geospatial professional services<br>
<a href="http://www.spatialys.com" target="_blank">http://www.spatialys.com</a><br>
_______________________________________________<br>
gdal-dev mailing list<br>
<a href="mailto:gdal-dev@lists.osgeo.org">gdal-dev@lists.osgeo.org</a><br>
<a href="http://lists.osgeo.org/mailman/listinfo/gdal-dev" target="_blank">http://lists.osgeo.org/mailman/listinfo/gdal-dev</a></div></div></blockquote></div><br></div>