On Thu, Jun 17, 2010 at 1:40 PM, Alexander Bruy <[email protected]> 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, [email protected] light and sound - activate the windows | http://pobox.com/~warmerdam and watch the world go round - Rush | Geospatial Programmer for Rent _______________________________________________ gdal-dev mailing list [email protected] http://lists.osgeo.org/mailman/listinfo/gdal-dev
