# [gdal-dev] Real and "raster" coordinates

Pinner, Luke Luke.Pinner at environment.gov.au
Wed Mar 24 18:34:46 EDT 2010

```The following module might help, use the mapToPixel() function.

#---Imports
try:
from osgeo import gdal
except ImportError:
import gdal

def mapToPixel(mx,my,gt):
''' Convert map to pixel coordinates
@param  mx    Input map x coordinate (double)
@param  my    Input map y coordinate (double)
@param  gt    Input geotransform (six doubles)
@return px,py Output coordinates (two doubles)
'''
if gt[2]+gt[4]==0: #Simple calc, no inversion required
px = (mx - gt[0]) / gt[1]
py = (my - gt[3]) / gt[5]
else:
px,py=ApplyGeoTransform(mx,my,InvGeoTransform(gt))
return int(px+0.5),int(py+0.5)

def pixelToMap(px,py,gt):
''' Convert pixel to map coordinates
@param  px    Input pixel x coordinate (double)
@param  py    Input pixel y coordinate (double)
@param  gt    Input geotransform (six doubles)
@return mx,my Output coordinates (two doubles)
'''
mx,my=ApplyGeoTransform(px,py,gt)
return mx,my

def ApplyGeoTransform(inx,iny,gt):
''' Apply a geotransform
@param  inx       Input x coordinate (double)
@param  iny       Input y coordinate (double)
@param  gt        Input geotransform (six doubles)
@return outx,outy Output coordinates (two doubles)
'''
outx = gt[0] + inx*gt[1] + iny*gt[2]
outy = gt[3] + inx*gt[4] + iny*gt[5]
return (outx,outy)

def InvGeoTransform(gt_in):
'''

************************************************************************
*                        InvGeoTransform(gt_in)

************************************************************************

**
* Invert Geotransform.
*
* This function will invert a standard 3x2 set of GeoTransform
coefficients.
*
* @param  gt_in  Input geotransform (six doubles - unaltered).
* @return gt_out Output geotransform (six doubles - updated) on
success,
*                None if the equation is uninvertable.
'''
#
************************************************************************
******
#    * This code ported from GDALInvGeoTransform() in
gdaltransformer.cpp
#    * as it isn't exposed in the python SWIG bindings until GDAL
1.7
#    * copyright & permission notices included below as per
conditions.
#
#
************************************************************************
******
#    * \$Id: gdaltransformer.cpp 15024 2008-07-24 19:25:06Z rouault \$
#    *
#    * Project:  Mapinfo Image Warper
#    * Purpose:  Implementation of one or more GDALTrasformerFunc
types, including
#    *           the GenImgProj (general image reprojector)
transformer.
#    * Author:   Frank Warmerdam, warmerdam at pobox.com
#    *
#
************************************************************************
******
#    * Copyright (c) 2002, i3 - information integration and imaging
#    *                          Fort Collin, CO
#    *
#    * Permission is hereby granted, free of charge, to any person
obtaining a
#    * copy of this software and associated documentation files (the
"Software"),
#    * to deal in the Software without restriction, including
without limitation
#    * the rights to use, copy, modify, merge, publish, distribute,
#    * and/or sell copies of the Software, and to permit persons to
whom the
#    * Software is furnished to do so, subject to the following
conditions:
#    *
#    * The above copyright notice and this permission notice shall
be included
#    * in all copies or substantial portions of the Software.
#    *
#    * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY
KIND, EXPRESS
#    * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
MERCHANTABILITY,
#    * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO
EVENT SHALL
#    * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,
DAMAGES OR OTHER
#    * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR
OTHERWISE, ARISING
#    * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
OTHER
#    * DEALINGS IN THE SOFTWARE.
#
************************************************************************
****

# we assume a 3rd row that is [1 0 0]

# Compute determinate
det = gt_in[1] * gt_in[5] - gt_in[2] * gt_in[4]

if( abs(det) < 0.000000000000001 ):
return

inv_det = 1.0 / det

# compute adjoint, and divide by determinate
gt_out = [0,0,0,0,0,0]
gt_out[1] =  gt_in[5] * inv_det
gt_out[4] = -gt_in[4] * inv_det

gt_out[2] = -gt_in[2] * inv_det
gt_out[5] =  gt_in[1] * inv_det

gt_out[0] = ( gt_in[2] * gt_in[3] - gt_in[0] * gt_in[5]) * inv_det
gt_out[3] = (-gt_in[1] * gt_in[3] + gt_in[0] * gt_in[4]) * inv_det

return gt_out

------
If you have received this transmission in error please notify us immediately by return e-mail and delete all copies. If this e-mail or any attachments have been sent to you in error, that error does not constitute waiver of any confidentiality, privilege or copyright in respect of information in the e-mail or attachments.

Please consider the environment before printing this email.

------

```