[gdal-dev] Need code review obtaining elev data in CPP Class (GDAL beginner)
Javier Jimenez Shaw
j1 at jimenezshaw.com
Sun Oct 24 15:58:55 PDT 2021
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 <lve0200 at gmail.com> 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
> gdal-dev at lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/gdal-dev
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20211025/168de017/attachment.html>
More information about the gdal-dev
mailing list