Latitude and Longitude to pixel-value in GDAL

Question

I'm trying to read a GeoTIFF with one band which stores a value between 0 and 3 ( 255 is no maks ). My Goal is to write a little program which takes in a latitude/longitude and returns the fitting pixel value at the geocoordinate in the geotiff.

I downloaded it from here : https://dataverse.harvard.edu/dataset.xhtml?persistentId=doi:10.7910/DVN/QQHCIK And if you wanna take a look at it... heres the download link. https://drive.google.com/file/d/104De_YQN1V8tSbj2uhIPO0FfowjnInnX/view?usp=sharing

However, my code does not work. I cant tell you what is wrong, it just outputs the wrong pixel value every damn time. I guess either my geocordinate to pixel transformation is wrong or theres a huge logic mistake.

This was my first try, it prints the wrong value for the geo-coordinate. Furthermore it crashes with some geocoordinates, a negative longitude makes it crashes because this will make pixelY negative which causes an exception in the raster method.

            GdalConfiguration.ConfigureGdal();
            var tif = "C:\\Users\\Lars\\Downloads\\pnv_biome.type_biome00k_c_1km_s0..0cm_2000..2017_v0.1.tif";
            var lat = 51.0;
            var lon = 8.3;
            using (var image = Gdal.Open(tif, Access.GA_ReadOnly)) {

                var bOneBand = image.GetRasterBand(1);
                var width = bOneBand.XSize;
                var height = bOneBand.YSize;
                
                // Geocoordinate ( lat, lon ) to pixel in the tiff ? 
                var geoTransform = new double[6];
                image.GetGeoTransform(geoTransform);
                Gdal.InvGeoTransform(geoTransform, geoTransform);
                
                var bOne = new int[1];
                Gdal.ApplyGeoTransform(geoTransform, lat, lon, out var pixelXF, out var pixelYF);

                var pixelX = (int)Math.Floor(pixelXF);
                var pixelY = (int)Math.Floor(pixelYF);
                
                // Read pixel 
                bOneBand.ReadRaster(pixelX, pixelY, 1, 1, bOne, 1, 1,0,0);
                Console.WriteLine(bOne[0]); // bOne[0] contains wrong data
            }

My second attempt looks like the following... but also outputs the wrong pixel value for the given coordinates. It also crashes with some geocoordinates.

GdalConfiguration.ConfigureGdal();
            var tif = "C:\\Users\\Lars\\Downloads\\pnv_biome.type_biome00k_c_1km_s0..0cm_2000..2017_v0.1.tif";
            var lat = 24.7377; // 131.5847
            var lon = 131.5847;
            
            using (var image = Gdal.Open(tif, Access.GA_ReadOnly)) {

                var bOneBand = image.GetRasterBand(1);
                var bOne = new int[1];
                
                // Spatial reference to transform latlng to map coordinates... from python code
                var point_srs = new SpatialReference("");
                var file_srs = new SpatialReference(image.GetProjection());
                point_srs.ImportFromEPSG(4326);
                point_srs.SetAxisMappingStrategy(AxisMappingStrategy.OAMS_TRADITIONAL_GIS_ORDER);
           
                var mapCoordinates = new double[2]{ lat, lon};
                var transform = new CoordinateTransformation(point_srs, file_srs);
                transform.TransformPoint(mapCoordinates);
                
                // Map coordinates to pixel coordinates ? 
                var geoTransform = new double[6];
                image.GetGeoTransform(geoTransform);
                Gdal.InvGeoTransform(geoTransform, geoTransform);
                
                Gdal.ApplyGeoTransform(geoTransform, mapCoordinates[0], mapCoordinates[1], out var pixelXF, out var pixelYF);
                var pixelX = (int)pixelXF;
                var pixelY = (int)pixelYF;
                
                bOneBand.ReadRaster(pixelX, pixelY, 1, 1, bOne, 1, 1, 0,0);
                Console.WriteLine(bOne[0]);  // bOne[0] contains wrong value
            }

What is wrong with my code, why does it output the wrong pixel value ? Any help appreciated !

Answer 1

Perhaps you are using lat/lon where it should be lon/lat? It would be helpful if you printed some values such as mapCoordinates.

Here is how you can do this with R, for comparison:

library(terra)
r = rast("pnv_biome.type_biome00k_c_1km_s0..0cm_2000..2017_v0.1.tif")
xy = cbind(c(8.3, 131.5847), c(51.0, 24.7377))
extract(r, xy)
#  pnv_biome.type_biome00k_c_1km_s0..0cm_2000..2017_v0.1
#1                                                     9
#2                                                    NA

And the zero-based row/col numbers

rowColFromCell(r, cellFromXY(r, xy)) - 1
#     [,1]  [,2]
#[1,] 4364 22596
#[2,] 7515 37390

And you can use describe (a.k.a. GDALinfo) to get some metadata. E.g.

d = describe("pnv_biome.type_biome00k_c_1km_s0..0cm_2000..2017_v0.1.tif")
d[50]
# [1] "Band 1 Block=43200x1 Type=Byte, ColorInterp=Gray"