14

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

4
  • 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, 2011 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, 2011 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, 2011 at 21:08
  • 2
    This (newer) thread provides an explicit example and another solution: gis.stackexchange.com/questions/8430/…
    – whuber
    Apr 12, 2011 at 15:35

3 Answers 3

11

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)

5
  • (+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, 2011 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. Apr 11, 2011 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, 2011 at 19:19
  • 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, 2011 at 19:47
  • According to the docs, the -a_srs parameter overrides the projection for the output file but it explicitly notes that no reprojection is done.
    – A. Bollans
    Apr 20, 2022 at 9:03
7

The geotransform is documented at https://gdal.org/user/raster_data_model.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.

4
  • 1
    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, 2011 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, 2011 at 19:45
  • Link is broken. @whuber am I right that the affine transform (from GDAL's GetGeoTransform) is a transform between pixels and CRS coordinates? So, for instance if the data is in a GDA projection, points need to be transformed from WGS84/lat/lon to GDAXX using e.g. the CoordianteTransform, and then to pixels using the affine GeoTransform? This two-step process is something that's been tripping me up quite a lot, and I suspect is causing the OP trouble too. Didn't help that xarray/rasterio appears to have a bug that screws up the GeoTransform values: github.com/pydata/xarray/issues/3185
    – naught101
    Jan 6, 2020 at 3:57
  • @naught I found the new URL and updated my post.
    – whuber
    Jan 6, 2020 at 15:47
0

You can use the following:

import gdal, numpy as n
def coord(file):
    padfTransform = file.GetGeoTransform()
    indices = n.indices(file.ReadAsArray().shape)
    xp = padfTransform[0] + indices[1]*padfTransform[1] + indices[1]*padfTransform[2]   
    yp = padfTransform[3] + indices[0]*padfTransform[4] + indices[0]*padfTransform[5]  
    return xp,yp
file = gdal.Open('Yourfile.tif') # A GeoTiff file
x,y = coord(file)

coord will return longitude(x) and latitude(y) of all the pixels. Keep in mind that the coordinates are of the left corner of a pixel

1
  • Shouldn't it be indixes[1]*padfTransform[1] + indices[0]*padfTransform[2] and vice versa for the yp? Apr 11, 2017 at 15:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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