2

I have written the following code that return the elevation given longitude and latitude from a GeoTIFF file:

    private static double getValue( String filePath, double longitude, double latitude) throws Exception {

        gdal.AllRegister();

        Dataset dataset = gdal.Open(filePath, gdalconstConstants.GA_ReadOnly);

        String projection = dataset.GetProjection();
        double[] transform = dataset.GetGeoTransform();

        SpatialReference src = new SpatialReference();
        src.SetWellKnownGeogCS("WGS84");
        SpatialReference dst = new SpatialReference(projection);
        CoordinateTransformation coordinateTransformation = new CoordinateTransformation(src, dst);


        double[] xy = coordinateTransformation.TransformPoint(longitude, latitude);

        float[] result = new float[1];
        int[] band = {1};
        
        int x = (int)(((xy[0] - transform[0]) / transform[1]));
        int y = (int)(((xy[1] - transform[3]) / transform[5]));
        int readRaster = dataset.ReadRaster(x, y, 1, 1, 1, 1, 6, result, band);
        if (readRaster == 3){
            throw new IndexOutOfBoundsException();
        }
        double elevation = result[0];
        return elevation;
    }


I am working with the following geoTiff file https://download.osgeo.org/geotiff/samples/gdal_eg/

and I compare the result with the gdallocationinfo command .

the problem is that I do not get the same x,y pixels

for example if I use:

"longitude":-117.631169, "latitude" : 33.722355

In the code I get x = 14 and y = 408 and elevation = 132

but if I use it with gdallocationinfo cea.tif -wgs84 -117.631169 33.722355 I get x = 16 and y = 408 and elevation = 90

Another issue that I had that I notice that the elevation of the lon lat define in the example are not the real elevation in an elevation finder tool https://www.daftlogic.com/sandbox-google-maps-find-altitude.htm which there it says 579.00 m

1
  • can somebody help please?
    – Samuel
    Commented Oct 27, 2021 at 8:02

1 Answer 1

1

I use the following code, based on https://gis.stackexchange.com/a/415337/219296, which should work just fine:

class ALOSTile():
    def __init__(self, tile_path):
        """
        Load a tile, and cache it. Individual tiles are relatively
        small, so for performance we preload the entire tile's
        elevation data into RAM.
        """
        self.tile_path = tile_path
        self.dataset = gdal.Open(tile_path, gdal.GA_ReadOnly)

        if self.dataset is None:
            raise Exception(f'Could not load GDAL file{tile_path}')

        self.grid = self.dataset.GetRasterBand(1).ReadAsArray()

    def lookup(self, lat, lon):
        """
        see https://gis.stackexchange.com/a/415337/219296
        """
        try:
            forward_transform = self.dataset.GetGeoTransform()
            reverse_transform = gdal.InvGeoTransform(forward_transform)
            x, y = [int(v) for v in gdal.ApplyGeoTransform(
                reverse_transform, lon, lat)]
            return int(self.grid[y][x])

        except:
            return None

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.