If I understand properly what you did... there is a missing transformation in valueAt
You want to query with lat-lon in wgs84, but the image is in ETRS89_ETRS_LAEA So you need to transform the lat-lon in wgs84 to x-y in ETRS89_ETRS_LAEA, and then to the pixel coordinates using the inverseTransform. El dom., 24 oct. 2021 18:00, Hans Kurscheidt <[email protected]> escribió: > 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 > > > _______________________________________________ > gdal-dev mailing list > [email protected] > https://lists.osgeo.org/mailman/listinfo/gdal-dev >
_______________________________________________ gdal-dev mailing list [email protected] https://lists.osgeo.org/mailman/listinfo/gdal-dev
