<div dir="auto"><div>If I understand properly what you did... there is a missing transformation in valueAt</div><div dir="auto"><br></div><div dir="auto">You want to query with lat-lon in wgs84, but the image is in <span style="font-family:sans-serif">ETRS89_ETRS_LAEA</span></div><div dir="auto"><font face="sans-serif">So you need to transform the lat-lon in wgs84 to x-y in </font><span style="font-family:sans-serif">ETRS89_ETRS_LAEA, and then to the pixel coordinates using the </span><span style="font-family:sans-serif">inverseTransform.</span><font face="sans-serif"><br></font><br><div class="gmail_quote" dir="auto"><div dir="ltr" class="gmail_attr">El dom., 24 oct. 2021 18:00, Hans Kurscheidt <<a href="mailto:lve0200@gmail.com" target="_blank" rel="noreferrer">lve0200@gmail.com</a>> escribió:<br></div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">I came to realize that I apparently sent HTML. I try now send text only.<br>
<br>
hk<br>
<br>
Am 24.10.2021 um 17:42 schrieb Hans Kurscheidt:<br>
><br>
> Dears, excuse in my very 1st post here, a bloody beginner's question <br>
> and probably behavior. I need every now and then an elevation value <br>
> from WGS84 coordinates _in my application_. What does<br>
><br>
> gdallocationinfo -valonly -wgs84 /home/hk/eu_dem_v11_E30N20.TIF -1.0 <br>
> 44.0 at console level.<br>
><br>
> I did read the source code and large parts of the GDAL documentation <br>
> and the best Class initialization I was able to come up with is this: <br>
> (I made few vars static to avoid compiler optimization)<br>
><br>
> class GDALRasterImage {<br>
> public:<br>
> GDALRasterImage(const char* filename);<br>
> ~GDALRasterImage();<br>
> int valueAt(double lat, double lon);<br>
><br>
> private:<br>
> GDALDatasetH dataset;<br>
> GDALRasterBandH band;<br>
> char *pszSourceSRS = nullptr;<br>
> double inverseTransform[6];<br>
> double datasetTransform[6];<br>
> };<br>
> GDALRasterImage::GDALRasterImage(const char* filename) {<br>
><br>
> double llon = -1, llat = 44; static bool rv;<br>
> GDALAllRegister();<br>
><br>
> // SanitizeSRS to WGS84<br>
> CPLErrorReset();<br>
> OGRSpatialReferenceH hSRS = OSRNewSpatialReference(nullptr);<br>
> CPLFree(pszSourceSRS);<br>
> if (OSRSetFromUserInput(hSRS, "WGS84") == OGRERR_NONE)<br>
> OSRExportToWkt(hSRS, &pszSourceSRS);<br>
> OSRDestroySpatialReference(hSRS);<br>
><br>
> this->dataset = GDALOpenEx(filename, GDAL_OF_RASTER, nullptr, <br>
> nullptr, nullptr);<br>
> assert(dataset != nullptr);<br>
><br>
> // Assume there is only one band in the raster source and use that<br>
> assert(GDALGetRasterCount(dataset) == 1);<br>
> this->band = GDALGetRasterBand(dataset, 1);<br>
><br>
> // Setup coordinate transformation<br>
> static OGRSpatialReferenceH hSrcSRS = nullptr, hTrgSRS = nullptr;<br>
> static OGRCoordinateTransformationH hCT = nullptr;<br>
> hSrcSRS = OSRNewSpatialReference(pszSourceSRS);<br>
> hTrgSRS = OSRNewSpatialReference(GDALGetProjectionRef(dataset));<br>
><br>
> hCT = OCTNewCoordinateTransformation(hSrcSRS, hTrgSRS);<br>
> if (hCT == nullptr) exit(1);<br>
> rv= OCTTransform(hCT, 1, &llon, &llat, nullptr);<br>
><br>
> // Get the inverse geo transform to map from geo location -> pixel <br>
> location<br>
> if (GDALGetGeoTransform(dataset, datasetTransform) != CE_None) <br>
> exit (1);<br>
> if (!GDALInvGeoTransform(datasetTransform, inverseTransform)) exit(1);<br>
> }<br>
> The init runs w/out obvious code failures, but when I call<br>
><br>
> int GDALRasterImage::valueAt(double lat, double lon) {<br>
> int x = static_cast<int>(floor(inverseTransform[0] + <br>
> inverseTransform[1] * lon + inverseTransform[2] * lat));<br>
> int y = static_cast<int>(floor(inverseTransform[3] + <br>
> inverseTransform[4] * lon + inverseTransform[5] * lat));<br>
><br>
> int32_t pixelValue;<br>
> GDALRasterIO(band, GF_Read, x, y, 1, 1, &pixelValue, 1, 1, <br>
> GDT_Int32, 0, 0);<br>
> return pixelValue;<br>
> }<br>
><br>
> with Lon=-1.0 and Lat=44.0 I receive a run time error:<br>
><br>
> ERROR 5: /home/hk/eu_dem_v11_E30N20.TIF, band 1: Access window out of <br>
> range in RasterIO(). Requested<br>
> (-120001,119998) of size 1x1 on raster of 40000x40000.<br>
><br>
> I guess it has to do with the inverse transformation, but for the life <br>
> of me, I have not the slightest idea, where to search from here.<br>
><br>
> gdalinfo /home/hk/eu_dem_v11_E30N20.TIF<br>
> Driver: GTiff/GeoTIFF<br>
> Files: /home/hk/eu_dem_v11_E30N20.TIF<br>
> Size is 40000, 40000<br>
> Coordinate System is:<br>
> PROJCS["<span style="font-family:sans-serif">ETRS89_ETRS_LAEA</span>",<br>
> GEOGCS["ETRS89",<br>
> DATUM["European_Terrestrial_Reference_System_1989",<br>
> SPHEROID["GRS 1980",6378137,298.2572221010042,<br>
> AUTHORITY["EPSG","7019"]],<br>
> AUTHORITY["EPSG","6258"]],<br>
> PRIMEM["Greenwich",0],<br>
> UNIT["degree",0.0174532925199433],<br>
> AUTHORITY["EPSG","4258"]],<br>
> PROJECTION["Lambert_Azimuthal_Equal_Area"],<br>
> PARAMETER["latitude_of_center",52],<br>
> PARAMETER["longitude_of_center",10],<br>
> PARAMETER["false_easting",4321000],<br>
> PARAMETER["false_northing",3210000],<br>
> UNIT["metre",1,<br>
> AUTHORITY["EPSG","9001"]]]<br>
> Origin = (3000000.000000000000000,3000000.000000000000000)<br>
> Pixel Size = (25.000000000000000,-25.000000000000000)<br>
> Metadata:<br>
> AREA_OR_POINT=Area<br>
> DataType=Elevation<br>
> Image Structure Metadata:<br>
> COMPRESSION=LZW<br>
> INTERLEAVE=BAND<br>
> Corner Coordinates:<br>
> Upper Left ( 3000000.000, 3000000.000) ( 8d 7'39.52"W, 48d38'23.47"N)<br>
> Lower Left ( 3000000.000, 2000000.000) ( 5d28'46.07"W, 39d52'33.70"N)<br>
> Upper Right ( 4000000.000, 3000000.000) ( 5d31' 4.07"E, 50d 1'26.83"N)<br>
> Lower Right ( 4000000.000, 2000000.000) ( 6d11'52.97"E, 41d 1'36.49"N)<br>
> Center ( 3500000.000, 2500000.000) ( 0d27' 3.13"W, 45d 5'35.85"N)<br>
> Band 1 Block=128x128 Type=Float32, ColorInterp=Gray<br>
> NoData Value=-3.4028234663852886e+38<br>
> Metadata:<br>
> BandName=Band_1<br>
> RepresentationType=ATHEMATIC<br>
><br>
> The transformation data I have are:<br>
><br>
> datasetTransform[6] contains: 3000000, 25, 0, 3000000, -25<br>
><br>
> inverseTransform[6] contains: -120000, 0.04, 0, 12000, -0.04, 0<br>
><br>
><br>
> Thanks for your help!<br>
><br>
> hans<br>
><br>
_______________________________________________<br>
gdal-dev mailing list<br>
<a href="mailto:gdal-dev@lists.osgeo.org" rel="noreferrer noreferrer" target="_blank">gdal-dev@lists.osgeo.org</a><br>
<a href="https://lists.osgeo.org/mailman/listinfo/gdal-dev" rel="noreferrer noreferrer noreferrer" target="_blank">https://lists.osgeo.org/mailman/listinfo/gdal-dev</a><br>
</blockquote></div></div></div>