[Gdal-dev] Get a pixel from a lat/lon.

[Frank Warmerdam]

Stephane Routelous wrote:
> hi,
> 
> funny, I'm just facing the same question !
> 
> what to you mean by 2) ?

   2) Apply the hTransform (with OCTTransform) before using the geotransform.

In the code Ken provided he was using the geotransform to convert from
pixel/line to projected coordinates, and then calling OCTTransform()
to reproject to lat/long.  When going from lat/long to pixel/line you
need to apply the OCTTransform() call to the input values and then do the
part with geotransform.

> And I'm also interessed by the 3) but with adfGeoTransform[2] and
> adfGeoTransform[4] != 0 . I tried to solve the equation by myself, my math
> courses are quite old :-(

The "properly" solution is to invert the geotransform and apply it the old way.
In a relatively recent GDAL the GDALInvGeoTransform() was made a public function.
So we would use instead:

     double adfInvGeoTransform[6];

     GDALInvGeoTransform( adfGeoTransform, adfInvGeoTransform );

     x = adfInvGeoTransform[0] + adfInvGeoTransform[1] * dfGeoX
                               + adfInvGeoTransform[2] * y;
     y = adfInvGeoTransform[3] + adfInvGeoTransform[4] * dfGeoY
                               + adfInvGeoTransform[5] * y;

The code for GDALInvGeoTransform() is in gdal/alg/gdaltransformer.cpp if you
want to inspect it.  It is a simple matrix inversion, but I have to admit my
math is also pretty rusty and I had to find the equations on the web.

You can also use GDALApplyGeoTransform() if you don't like the explicit use of
the geotransform equations by the way (make sure x and y are doubles!):

     GDALApplyGeoTransform( adfInvGeoTransform, dfGeoX, dfGeoY, &x, &y );

Best regards,

-- 
---------------------------------------+--------------------------------------
I set the clouds in motion - turn up   | Frank Warmerdam, warmerdam at pobox.com
light and sound - activate the windows | http://pobox.com/~warmerdam
and watch the world go round - Rush    | Geospatial Programmer for Rent