[gdal-dev] Real and "raster" coordinates

[Pinner, Luke]

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,
sublicense,
    #    * 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.

------