[Gdal-dev] UTM to Lat/Lon

[Andrey Kiselev]

On Thu, Jul 03, 2003 at 11:10:14AM +0200, mathieu_gdal wrote:
> I would like to implement to raster projection: Converting a pixel
> coordinate (in the bitmap) into a position in Lat/Lon coord system (in
> WGS84 for example)
> 
> I guess the way to do it is:
> 
> 1) Convert from pixel coordinate to the geotiff file coordinate system,
> ie in my case UTM zone 19, in meters. (There is no problem for that, I
> will use the affine coefficients)
> 
> 2) Convert from the geotiff file coordinate system to lat/lon system.
> 
> (u,v)Pixel ---------1----------> UTM(x,y) meters ---------2---------->
> (Lat,Lon)radians
> 
> How to do 2) ??

There is a simple Python script attached to read geographic information
from the GDAL supported file and convert specified pixel/line
coordinates into lat/long ones. It is pretty self-explanatory and will
be added to gdal/pymod/samples soon.

> Note : I have tried to create a OGRCoordinateTransformation with my
> file's OGRSpatialReference and a simple lat/lon OGRSpatialReference. It
> does not work, even if the transformation seems to be well created.

It works ;-) If you need working C code, look into gdailnfo.c sources,
GDALInfoReportCorner() function.


-- 
Andrey V. Kiselev
Home phone:  +7 812 5274898  ICQ# 26871517
-------------- next part --------------
#!/usr/bin/env python

import gdal
from gdalconst import *
import osr
import sys

# =============================================================================
def Usage():
    print
    print 'Read coordinate system and geotransformation matrix from input'
    print 'file and print latitude/longitude coordinates for the specified'
    print 'pixel.'
    print
    print 'Usage: tolatlong.py pixel line infile'
    print
    sys.exit( 1 )

# =============================================================================

infile = None
pixel = None
line = None

# =============================================================================
# Parse command line arguments.
# =============================================================================
i = 1
while i < len(sys.argv):
    arg = sys.argv[i]

    if pixel is None:
	pixel = float(arg)

    elif line is None:
	line = float(arg)

    elif infile is None:
	infile = arg

    else:
	Usage()

    i = i + 1

if infile is None:
    Usage()
if pixel is None:
    Usage()
if line is None:
    Usage()

indataset = gdal.Open( infile, GA_ReadOnly )

geomatrix = indataset.GetGeoTransform()
X = geomatrix[0] + geomatrix[1] * pixel + geomatrix[2] * line
Y = geomatrix[3] + geomatrix[4] * pixel + geomatrix[5] * line

srs = osr.SpatialReference()
srs.ImportFromWkt(indataset.GetProjection())

srsLatLong = srs.CloneGeogCS()
ct = osr.CoordinateTransformation(srs, srsLatLong)
(long, lat, height) = ct.TransformPoint(X, Y)

print 'pixel: %g, line: %g' % (pixel, line)
print 'latitude: %g, longitude: %g (in degrees)' % (long, lat)