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 _in my application_. What does
gdallocationinfo -valonly -wgs84 /home/hk/eu_dem_v11_E30N20.TIF -1.0
44.0 at console level.
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)
class GDALRasterImage {
public:
GDALRasterImage(const char* filename);
~GDALRasterImage();
int valueAt(double lat, double lon);
private:
GDALDatasetH dataset;
GDALRasterBandH band;
char *pszSourceSRS = nullptr;
double inverseTransform[6];
double datasetTransform[6];
};
GDALRasterImage::GDALRasterImage(const char* filename) {
double llon = -1, llat = 44; static bool rv;
GDALAllRegister();
// SanitizeSRS to WGS84
CPLErrorReset();
OGRSpatialReferenceH hSRS = OSRNewSpatialReference(nullptr);
CPLFree(pszSourceSRS);
if (OSRSetFromUserInput(hSRS, "WGS84") == OGRERR_NONE)
OSRExportToWkt(hSRS, &pszSourceSRS);
OSRDestroySpatialReference(hSRS);
this->dataset = GDALOpenEx(filename, GDAL_OF_RASTER, nullptr,
nullptr, nullptr);
assert(dataset != nullptr);
// Assume there is only one band in the raster source and use that
assert(GDALGetRasterCount(dataset) == 1);
this->band = GDALGetRasterBand(dataset, 1);
// Setup coordinate transformation
static OGRSpatialReferenceH hSrcSRS = nullptr, hTrgSRS = nullptr;
static OGRCoordinateTransformationH hCT = nullptr;
hSrcSRS = OSRNewSpatialReference(pszSourceSRS);
hTrgSRS = OSRNewSpatialReference(GDALGetProjectionRef(dataset));
hCT = OCTNewCoordinateTransformation(hSrcSRS, hTrgSRS);
if (hCT == nullptr) exit(1);
rv= OCTTransform(hCT, 1, &llon, &llat, nullptr);
// Get the inverse geo transform to map from geo location -> pixel
location
if (GDALGetGeoTransform(dataset, datasetTransform) != CE_None)
exit (1);
if (!GDALInvGeoTransform(datasetTransform, inverseTransform)) exit(1);
}
The init runs w/out obvious code failures, but when I call
int GDALRasterImage::valueAt(double lat, double lon) {
int x = static_cast<int>(floor(inverseTransform[0] +
inverseTransform[1] * lon + inverseTransform[2] * lat));
int y = static_cast<int>(floor(inverseTransform[3] +
inverseTransform[4] * lon + inverseTransform[5] * lat));
int32_t pixelValue;
GDALRasterIO(band, GF_Read, x, y, 1, 1, &pixelValue, 1, 1,
GDT_Int32, 0, 0);
return pixelValue;
}
with Lon=-1.0 and Lat=44.0 I receive a run time error:
ERROR 5: /home/hk/eu_dem_v11_E30N20.TIF, band 1: Access window out of
range in RasterIO(). Requested
(-120001,119998) of size 1x1 on raster of 40000x40000.
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.
gdalinfo /home/hk/eu_dem_v11_E30N20.TIF
Driver: GTiff/GeoTIFF
Files: /home/hk/eu_dem_v11_E30N20.TIF
Size is 40000, 40000
Coordinate System is:
PROJCS["ETRS89_ETRS_LAEA",
GEOGCS["ETRS89",
DATUM["European_Terrestrial_Reference_System_1989",
SPHEROID["GRS 1980",6378137,298.2572221010042,
AUTHORITY["EPSG","7019"]],
AUTHORITY["EPSG","6258"]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433],
AUTHORITY["EPSG","4258"]],
PROJECTION["Lambert_Azimuthal_Equal_Area"],
PARAMETER["latitude_of_center",52],
PARAMETER["longitude_of_center",10],
PARAMETER["false_easting",4321000],
PARAMETER["false_northing",3210000],
UNIT["metre",1,
AUTHORITY["EPSG","9001"]]]
Origin = (3000000.000000000000000,3000000.000000000000000)
Pixel Size = (25.000000000000000,-25.000000000000000)
Metadata:
AREA_OR_POINT=Area
DataType=Elevation
Image Structure Metadata:
COMPRESSION=LZW
INTERLEAVE=BAND
Corner Coordinates:
Upper Left ( 3000000.000, 3000000.000) ( 8d 7'39.52"W, 48d38'23.47"N)
Lower Left ( 3000000.000, 2000000.000) ( 5d28'46.07"W, 39d52'33.70"N)
Upper Right ( 4000000.000, 3000000.000) ( 5d31' 4.07"E, 50d 1'26.83"N)
Lower Right ( 4000000.000, 2000000.000) ( 6d11'52.97"E, 41d 1'36.49"N)
Center ( 3500000.000, 2500000.000) ( 0d27' 3.13"W, 45d 5'35.85"N)
Band 1 Block=128x128 Type=Float32, ColorInterp=Gray
NoData Value=-3.4028234663852886e+38
Metadata:
BandName=Band_1
RepresentationType=ATHEMATIC
The transformation data I have are:
datasetTransform[6] contains: 3000000, 25, 0, 3000000, -25
inverseTransform[6] contains: -120000, 0.04, 0, 12000, -0.04, 0
Thanks for your help!
hans
