[gdal-dev] Problem with pixel/line coordinates

[Frank Warmerdam]

On Thu, Jun 17, 2010 at 1:40 PM, Alexander Bruy
<alexander.bruy at gmail.com> wrote:
> Hi,
>
> I use next code to translate real coordinates (lat/lon) to the pixel/line
> coordinates
>
> def mapToPixel( mX, mY, geoTransform ):
>  if geoTransform[ 2 ] + geoTransform[ 4 ] == 0:
>    pX = ( mX - geoTransform[ 0 ] ) / geoTransform[ 1 ]
>    pY = ( mY - geoTransform[ 3 ] ) / geoTransform[ 5 ]
>  else:
>    pX, pY = applyGeoTransform( mX, mY, invertGeoTransform( geoTransform ) )
>  return int( pX + 0.5 ), int( pY + 0.5 )

Alexander,

I believe the addition of 0.5 in the above is incorrect.  In the
simple, non-rotated case, all values from geoTransform[0]
to geoTransform[0] + geoTransform[1] would be on the 1st
pixel (ie. pX = 0).  As you have coded it, when you are half
way across the pixel you are switching into the next one.

Keep in mind that (geoTransform[0],geoTransform[3])
is the top left corner of the top left pixel - not the center.

> def pixelToMap( pX, pY, geoTransform ):
>  mX, mY = applyGeoTransform( pX, pY, geoTransform )
>  return mX, mY

Conversely, here if pX and pY are coming in as integer
rather than floating point values, then you will likely want
to add half a pixel before transforming so that you get
the geoeferenced location of the center of the pixel rather
than the upper left corner.

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