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

[Ken Sewell]

Frank,
	Thanks for your help, the GDALInvGeoTransform and
GDALApplyGeoTransform worked.  One quick correction, shouldn't the formula
below read like this:

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



Ken Sewell                   Defense Research Associates
Project Engineer             3915 Germany Lane, Suite 102
937-656-4899 x3004           Beavercreek, Ohio 45431-1608
sewell at dramail.com

On Thu, 30 Sep 2004, Frank Warmerdam wrote:

> 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
>
>