<html>
<head>
<meta http-equiv="content-type" content="text/html; charset=UTF-8">
</head>
<body>
<p>Dears, excuse in my very 1st post here, a bloody beginner's
question and probably behavior. I need every now and then an
elevation value from WGS84 coordinates <u>in my application</u>.
What does <br>
</p>
<p>gdallocationinfo -valonly -wgs84 /home/hk/eu_dem_v11_E30N20.TIF
-1.0 44.0 at console level.<br>
</p>
<p>I did read the source code and large parts of the GDAL
documentation and the best Class initialization I was able to come
up with is this: (I made few vars static to avoid compiler
optimization)<br>
</p>
<p>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, 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 location<br>
if (GDALGetGeoTransform(dataset, datasetTransform) != CE_None)
exit (1);<br>
if (!GDALInvGeoTransform(datasetTransform, inverseTransform))
exit(1);<br>
}<br>
The init runs w/out obvious code failures, but when I call</p>
<p>int GDALRasterImage::valueAt(double lat, double lon) {<br>
int x = static_cast<int>(floor(inverseTransform[0] +
inverseTransform[1] * lon + inverseTransform[2] * lat));<br>
int y = static_cast<int>(floor(inverseTransform[3] +
inverseTransform[4] * lon + inverseTransform[5] * lat));<br>
<br>
int32_t pixelValue;<br>
GDALRasterIO(band, GF_Read, x, y, 1, 1, &pixelValue, 1, 1,
GDT_Int32, 0, 0);<br>
return pixelValue;<br>
}</p>
<p>with Lon=-1.0 and Lat=44.0 I receive a run time error:</p>
<p>ERROR 5: /home/hk/eu_dem_v11_E30N20.TIF, band 1: Access window
out of range in RasterIO(). Requested<br>
(-120001,119998) of size 1x1 on raster of 40000x40000.<br>
</p>
<p>I guess it has to do with the inverse transformation, but for the
life of me, I have not the slightest idea, where to search from
here.</p>
<p>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["ETRS89_ETRS_LAEA",<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</p>
<p>The transformation data I have are:<br>
</p>
<p>datasetTransform[6] contains: 3000000, 25, 0, 3000000, -25</p>
<p>inverseTransform[6] contains: -120000, 0.04, 0, 12000, -0.04, 0</p>
<p><br>
</p>
<p>Thanks for your help!</p>
<p>hans<br>
</p>
</body>
</html>