[gdal-dev] Need code review obtaining elev data in CPP Class (GDAL beginner)

[Hans Kurscheidt]

I came to realize that I apparently sent HTML. I try now send text only.

hk

Am 24.10.2021 um 17:42 schrieb Hans Kurscheidt:
>
> 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
>