# How do I convert affine coordinates to lat/lng?

I am extremely new to GIS.

I'm using `gdal` to read in a landuse/landcover map and I need to pick out the lat/lng of certain land cover types to index into a different dataset which is expressed only in lat/lng. Unfortuantely, I don't understand the form of the x and y coordinates given to me from the geotransform, specifically the `originX` and `originY` below:

``````geotransform = dataset.GetGeoTransform()
originX = geotransform[0]
originY = geotransform[3]
``````

Printing these values gives me coordinates like `(447466.693808, 4952570.40529)`. How do these relate to the original latitude and longitude?

Edit:

Here's a simple python example that got me what I was looking for:

``````srs = osr.SpatialReference()
srs.ImportFromWkt(dataset.GetProjection())

srsLatLong = srs.CloneGeogCS()
ct = osr.CoordinateTransformation(srs,srsLatLong)
print ct.TransformPoint(originX,originY)
``````

Stolen from: tolatlong.py

• It looks like your data are projected (eg. UTM) and you will need to know what the projection is in order to "unproject" it back to long/lat coordinates. – user681 Apr 11 '11 at 18:35
• @Dan Thanks, so I know I can get the projection via a `dataset.GetProjectionRef()` and find out I'm using "UTM Zone 10", but then what? I'm googling around for methods like "unproject" but am coming up null. – Rich Apr 11 '11 at 18:41
• sorry for the use of the term unproject (in quotes) since data in decimal degrees is not projected, but if you want to get projected data back to decimal degrees from any given projection, then you have to (note the quotes) "project" it back to a geographic coordinate system, aka decimal degree data. – user681 Apr 11 '11 at 21:08
• This (newer) thread provides an explicit example and another solution: gis.stackexchange.com/questions/8430/… – whuber Apr 12 '11 at 15:35

gdal_translate will reproject your data from whatever projection it is in to any other (in this case you want EPSG:4326) using:

``````gdal_translate -a_srs epsg:4326 srcfile new_file
``````

or you could use gdaltrasform to convert the points (and I'm sure you can access that from Python(?) too)

• (+1) I don't know why this was downvoted, Ian, because it looks to me like some sort of operation like it is going to be needed after the affine transformation is applied. – whuber Apr 11 '11 at 19:07
• gdal_translate won't work the way you have it (or rather it would merely assign EPSG:4326 to the data), you need to warp the data with something like: gdalwarp -s_srs EPSG:32610 -t_srs EPSG:4326 src dest If the data has proj metadata already, you can ditch the -s_srs parameter. – MerseyViking Apr 11 '11 at 19:08
• Maybe this is what I'm after. Is the "epsg:4326" a transform to global lat/lng? I think I see some classes where I can provide a projection, I was thinking "WGS:84", but I don't know the difference. – Rich Apr 11 '11 at 19:19
• @Rich Click on the "coordinate-system" tag in your question, sort the resulting page by votes, and start reading. Many of your questions will be answered within a few minutes. – whuber Apr 11 '11 at 19:47

The geotransform is documented at http://www.gdal.org/gdal_datamodel.html . The idea is that you take (x,y) coordinates straight from the dataset, apply a linear transformation to get (u,v) with

``````u = a*x + b*y
v = c*x + d*y
``````

(you can take this to be the definition of a linear transformation), then shift the origin by adding geotransform[0] to u and geotransform[3] to v. That gives the "affine transformation" of (x,y). It's really intended to rotate, change scale, perhaps correct a little for some skew errors, and re-situate the data-specific coordinates (x,y) to match a known coordinate system. The result is supposed to produce projected coordinates. This simply means there is a mathematical procedure taking (longitude, latitude) and turning them into the known coordinates: this is called the "projection." "Unprojecting" is doing the reverse; so, if you know which projection is needed, you apply that to the affine transformed (x,y) coordinates to get the latitude and longitude.

By the way, the values of the constants a, b, c, d are given by the 1, 2, 4, and 5 entries in the geotransform array.

• I read the geotransform section of the link you provided, but I don't think this is what I'm after. Sorry if I am being obtuse, but doesn't this describe how to project x and y indices (0 to XSize or YSize) to projected coordinates? I'm curious how you translate projected coordinates into latitude and longitude. As a simple example, the geotransform defines the x/y origin of the upper left hand corner of the raster in projected coordinates (447466.693808, 4952570.40529). How would I turn those coordinates back into lat/lng (something like lat: 44,lng: -122). – Rich Apr 11 '11 at 19:17
• @Rich No, that link does not describe a projection. It only describes a change of coordinates between two maps, not between the world and a map. The unprojection turns the affine transform of the coordinates into (lat, lon). You do not want to write the code for that if you can help it! Unprojection is almost always a matter of identifying what projection is needed and calling the right software to do the job for you (as suggested in @iant's reply). – whuber Apr 11 '11 at 19:45

You can use the following:

``````import gdal, numpy as n
def coord(file):