[GIS] Convert latitude longitude pair to pixels in geotiff

coordinate systemcoordinatesgdal

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?

Best Answer

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.

Related Question