Convert latitude longitude pair to pixels in geotiff

Question

I am not a GIS guy.

I have a GeoTiff and I need to read elevation information out of it using GDAL (in some Java Code). So now I need to convert a given latitude/longitude pair into a pixel inside the GeoTiff.

Here's some information gdalinfo gave me about the file:

Driver: GTiff/GeoTIFF
Size is 58808, 30323
Coordinate System is:
PROJCS["MGI_Austria_Lambert",
    GEOGCS["GCS_MGI",
        DATUM["Militar_Geographische_Institute",
            SPHEROID["Bessel_1841",6377397.155,299.1528128,
                AUTHORITY["EPSG","7004"]],
            AUTHORITY["EPSG","6312"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Lambert_Conformal_Conic_2SP"],
    PARAMETER["standard_parallel_1",46],
    PARAMETER["standard_parallel_2",49],
    PARAMETER["latitude_of_origin",47.5],
    PARAMETER["central_meridian",13.33333333333333],
    PARAMETER["false_easting",400000],
    PARAMETER["false_northing",400000],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]
Origin = (106549.267203768889885,576922.512073625810444)
Pixel Size = (10.000000000000000,-10.000000000000000)
Metadata:
AREA_OR_POINT=Area
Image Structure Metadata:
COMPRESSION=LZW
INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  106549.267,  576922.512) (  9d19' 7.63"E, 49d 1'24.32"N)
Lower Left  (  106549.267,  273692.512) (  9d31'19.50"E, 46d17'55.02"N)
Upper Right (  694629.267,  576922.512) ( 17d21'50.31"E, 49d 1'22.34"N)
Lower Right (  694629.267,  273692.512) ( 17d 9'35.51"E, 46d17'53.15"N)
Center      (  400589.267,  425307.512) ( 13d20'28.29"E, 47d43'39.80"N)
Band 1 Block=128x128 Type=Float32, ColorInterp=Gray
NoData Value=-3.40282306073709653e+38

Can somebody tell me the formula on how to convert my coordinates into pixels?

Answer 1

Turns out that gdal provides everything needed (which should be no big surprise as gdal can do almost EVERYTHING :) )

So I just need to use the SpatialReference class to create my source projection which is WGS84:

SpatialReference src = new SpatialReference();
src.SetWellKnownGeogCS("WGS84");

and then my target projection:

dataset = gdal.Open("path/to/my/file", gdalconstConstants.GA_ReadOnly);
projection = dataset.GetProjection();
SpatialReference dst = new SpatialReference(projection);

Now I create a Transformation object like so:

CoordinateTransformation ct = new CoordinateTransformation(src, dst);

And then I transform a given pair of lat/lon to a pixel in my geotiff like so:

double[] xy = ct.TransformPoint(lon, lat);

int x = (int)(((xy[0] - transform[0]) / transform[1]));
int y = (int)(((xy[1] - transform[3]) / transform[5]));

You could further improve my taking a mean value of surrounding pixels as well, but for me that's enough.