Re: [gdal-dev] Problem with pixel/line coordinates
Thanks for your explanation, Frank! Now it works fine On Thu, 17 Jun 2010 13:59:46 -0400 Frank Warmerdam warmer...@pobox.com wrote: 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. -- Alexander Bruy ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] Problem with pixel/line coordinates
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.001: 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 ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Problem with pixel/line coordinates
On Thu, Jun 17, 2010 at 1:40 PM, Alexander Bruy alexander.b...@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, warmer...@pobox.com 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 gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
