1

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?

  • First create a point with the lat/lon, set the geographic spatial reference then use transformTo gdal.org/classOGRGeometry.html to convert the point into the spatial reference then find the offset in cells from the GeoTransform gdal.org/… then with RasterIO read the 1 by 1 array from the raster. I'm not a Java programmer but have done this in python, C# and C++ pretty much the same way with only a few minor changes for each language. – Michael Stimson Jun 27 '17 at 5:45
  • Thanks for this great hint! I think I got it working using the SpatialReference as well as the CoordinateTransformation class! Thanks! Do you want wo copy your comment as an answer so that I can accept it? – Georg Jun 27 '17 at 6:16
  • It would be better if you answered your own question, like I said I'm not a Java programmer so would struggle to put in a working code block. There is a meta post gis.meta.stackexchange.com/questions/4608/… decrying the answering of questions with a different software package so if I answered in C++, C# or python it would likely receive downvotes - and rightly so. – Michael Stimson Jun 27 '17 at 20:51
4

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.

  • +1 from me. You can read the pixels as an array and calculate the mean very quickly, the Java I've seen uses ReadRaster gdal.org/java/org/gdal/gdal/Band.html but be careful not to read rows or columns outside the pixel bounds. That's worthy of a new question though. – Michael Stimson Jun 28 '17 at 20:50
  • int x = (int)(((xy[0] - transform[0]) / transform[1])); What is the transform array here? I don't see any such variable in the code that you have posted. – Nancy Apr 1 at 9:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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