[Gdal-dev] How to convert Lat/Long -> Pixel/Line ?

[Frank Warmerdam]

Jérôme MAUPAY wrote:
> Sorry for that simple following question, but I really don't find the 
> answer in the doc:
>  
> I convert a pixel/line to lat/long using the following code:
>  
> bool CoxGisGdalFile::PixelLineToLatLong(int aPixel,int aLine,double* 
> aLat, double* aLong)
> {
>  // Cf http://www.remotesensing.org/gdal/class_GDALDataset.html#a8
>   *aLat = mGeoTransform[3] + mGeoTransform[4] * (double)aPixel + 
> mGeoTransform[5] * (double)aLine ;
>   *aLong = mGeoTransform[0] + mGeoTransform[1] * (double)aPixel + 
> mGeoTransform[2] * (double)aLine ;
>  
>     if( !OCTTransform(mTransform,1,aLong,aLat,NULL) )
>         return false ;
>     return true ;
> }
> But how to convert from lat/long to pixel/line ???

Jerome,

I realize this isn't a very fast answer.  I use the following function
to "invert" the geotransform when I want to translate the other way.

/************************************************************************/
/*                          InvGeoTransform()                           */
/*                                                                      */
/*      Invert a standard 3x2 "GeoTransform" style matrix with an       */
/*      implicit [1 0 0] final row.                                     */
/************************************************************************/

static int InvGeoTransform( double *gt_in, double *gt_out )

{
     double	det, inv_det;

     /* we assume a 3rd row that is [1 0 0] */

     /* Compute determinate */

     det = gt_in[1] * gt_in[5] - gt_in[2] * gt_in[4];

     if( fabs(det) < 0.000000000000001 )
         return 0;

     inv_det = 1.0 / det;

     /* compute adjoint, and devide by determinate */

     gt_out[1] =  gt_in[5] * inv_det;
     gt_out[4] = -gt_in[4] * inv_det;

     gt_out[2] = -gt_in[2] * inv_det;
     gt_out[5] =  gt_in[1] * inv_det;

     gt_out[0] = ( gt_in[2] * gt_in[3] - gt_in[0] * gt_in[5]) * inv_det;
     gt_out[3] = (-gt_in[1] * gt_in[3] + gt_in[0] * gt_in[4]) * inv_det;

     return 1;
}


You can see this used in GDALCreateGenImgProjTransformer() function in
gdal/alg/gdaltransformer.cpp (in CVS at least) which is documented at:

   http://www.remotesensing.org/gdal/gdal_alg_h.html#a4

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