[gdal-dev] Problem with pixel/line coordinates
Frank Warmerdam
warmerdam at pobox.com
Thu Jun 17 13:59:46 EDT 2010
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
More information about the gdal-dev
mailing list