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

Samuel
  • 21
  • 2

1 Answers1

0

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
Mike 'Pomax' Kamermans
  • 133
  • 7