[gdal-dev] Real and "raster" coordinates

classic Classic list List threaded Threaded
5 messages Options
Reply | Threaded
Open this post in threaded view
|

[gdal-dev] Real and "raster" coordinates

Alexander Bruy
Hi,

I have a georeferenced raster and coordinates of some point (lat, lon)
inside raster. Is it possible
to convert this lat/lon coordinates to raster coordinates (column and row)?

I use GDAL 1.6.3 with Python 2.5 (installed with OSGeo4W installer).

--
Alexander Bruy
_______________________________________________
gdal-dev mailing list
[hidden email]
http://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: Real and "raster" coordinates

chaitanya_ch
Alexander,

You can use the GDALGetGeoTransform() or GDALDataset::GetGeoTransform() to get the affine transformation coefficients which can be used to calculate lat/long to col/row and vice versa. You can see the utility script gdal_merge.py for an example.

On Wed, Mar 24, 2010 at 11:57 AM, Alexander Bruy <[hidden email]> wrote:
Hi,

I have a georeferenced raster and coordinates of some point (lat, lon)
inside raster. Is it possible
to convert this lat/lon coordinates to raster coordinates (column and row)?

I use GDAL 1.6.3 with Python 2.5 (installed with OSGeo4W installer).

--
Alexander Bruy
_______________________________________________
gdal-dev mailing list
[hidden email]
http://lists.osgeo.org/mailman/listinfo/gdal-dev



--
Best regards,
Chaitanya kumar CH.
/tʃaɪθənjə/ /kʊmɑr/
+91-9848167848
17.2416N 80.1426E

_______________________________________________
gdal-dev mailing list
[hidden email]
http://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: Real and "raster" coordinates

Jose Gomez-Dans
In reply to this post by Alexander Bruy
Hi

On 24 March 2010 06:27, Alexander Bruy <[hidden email]> wrote:
I have a georeferenced raster and coordinates of some point (lat, lon)
inside raster. Is it possible
to convert this lat/lon coordinates to raster coordinates (column and row)?

Use the geotransform (you may need to reproject your lat/lon into whatever projection your raster file is).

J

_______________________________________________
gdal-dev mailing list
[hidden email]
http://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

RE: Real and "raster" coordinates

lpinner
In reply to this post by Alexander Bruy
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, [hidden email]
    #    *
    #
************************************************************************
******
    #    * 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.

------

_______________________________________________
gdal-dev mailing list
[hidden email]
http://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: Real and "raster" coordinates

Alexander Bruy
Thanks! That what I need

2010/3/25 Pinner, Luke <[hidden email]>:

> 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):

--
Alexander Bruy
_______________________________________________
gdal-dev mailing list
[hidden email]
http://lists.osgeo.org/mailman/listinfo/gdal-dev