Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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

share|improve this question
    
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. –  Dan Patterson 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. –  Dan Patterson Apr 11 '11 at 21:08
1  
This (newer) thread provides an explicit example and another solution: gis.stackexchange.com/questions/8430/… –  whuber Apr 12 '11 at 15:35
add comment

2 Answers 2

up vote 5 down vote accepted

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)

share|improve this answer
    
(+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
    
Thanks iant, this put me on the right track. –  Rich Apr 11 '11 at 19:44
1  
@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
show 1 more comment

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.

share|improve this answer
    
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
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.