[gdal-dev] Problem with pixel/line coordinates

[Alexander Bruy]

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 )

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

def applyGeoTransform( inX, inY, geoTransform ):
  outX = geoTransform[ 0 ] + inX * geoTransform[ 1 ] + inY * geoTransform[ 2 ]
  outY = geoTransform[ 3 ] + inX * geoTransform[ 4 ] + inY * geoTransform[ 5 ]
  return outX, outY

def invertGeoTransform( geoTransform ):
  det = geoTransform[ 1 ] * geoTransform[ 5 ] - geoTransform[ 2 ] * geoTransform[ 4 ]

  if abs( det ) < 0.000000000000001:
    return

  invDet = 1.0 / det

  # compute adjoint and divide by determinate
  outGeoTransform = [ 0, 0, 0, 0, 0, 0 ]
  outGeoTransform[ 1 ] = geoTransform[ 5 ] * invDet
  outGeoTransform[ 4 ] = -geoTransform[ 4 ] * invDet

  outGeoTransform[ 2 ] = -geoTransform[ 2 ] * invDet
  outGeoTransfrom[ 5 ] = geoTransform[ 1 ] * invDet

  outGeoTransform[ 0 ] = ( geoTransform[ 2 ] * geoTransform[ 3 ] - geoTransform[ 0 ] * geoTransform[ 5 ] ) * invDet
  outGeoTransform[ 3 ] = ( -geoTransform[ 1 ] * geoTransform[ 3 ] + geoTransform[ 0 ] * geoTransform[ 4 ] ) * invDet

  return outGeoTransform

When I test it on large raster (in AIG/Arc/Info Binary Grid format)
I get next problem: raster value in some points reported by script
is differ from value reported by Info tool in QGIS and/or ArcGIS.
But for another points reported raster value is correct.

As I understand, when point is near pixel border, MapToPixel function
return wrong result - row and column of neighbour pixel. For example,
if point is near right side of the pixel I get column and row for right neighbour pixel not for pixel I need. Same problem when point is near
top side of pixel, I get colunm and row for top neighbour pixel.

Is it possible to get correct result for all points? May be there is a
mistake in code. Any help and suggestions are welcome.

I can upload sample data if necessary but they are large (~330 Mb)

I'm use GDAL 1.6.3 with Python 2.5.2 under Linux.

Thanks!

-- 
Alexander Bruy