python - Obtain Latitude and Longitude from a GeoTIFF File

Question

Using GDAL in Python, how do you get the latitude and longitude of a GeoTIFF file?

GeoTIFF's do not appear to store any coordinate information. Instead, they store the XY Origin coordinates. However, the XY coordinates do not provide the latitude and longitude of the top left corner and bottom left corner.

It appears I will need to do some math to solve this problem, but I don't have a clue on where to start.

What procedure is required to have this performed?

I know that the GetGeoTransform() method is important for this, however, I don't know what to do with it from there.

Answer 1

To get the coordinates of the corners of your geotiff do the following:

from osgeo import gdal
ds = gdal.Open('path/to/file')
width = ds.RasterXSize
height = ds.RasterYSize
gt = ds.GetGeoTransform()
minx = gt[0]
miny = gt[3] + width*gt[4] + height*gt[5] 
maxx = gt[0] + width*gt[1] + height*gt[2]
maxy = gt[3] 

However, these might not be in latitude/longitude format. As Justin noted, your geotiff will be stored with some kind of coordinate system. If you don't know what coordinate system it is, you can find out by running gdalinfo:

gdalinfo ~/somedir/somefile.tif 

Which outputs:

Driver: GTiff/GeoTIFF
Size is 512, 512
Coordinate System is:
PROJCS["NAD27 / UTM zone 11N",
    GEOGCS["NAD27",
        DATUM["North_American_Datum_1927",
            SPHEROID["Clarke 1866",6378206.4,294.978698213901]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",-117],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1]]
Origin = (440720.000000,3751320.000000)
Pixel Size = (60.000000,-60.000000)
Corner Coordinates:
Upper Left  (  440720.000, 3751320.000) (117d38'28.21"W, 33d54'8.47"N)
Lower Left  (  440720.000, 3720600.000) (117d38'20.79"W, 33d37'31.04"N)
Upper Right (  471440.000, 3751320.000) (117d18'32.07"W, 33d54'13.08"N)
Lower Right (  471440.000, 3720600.000) (117d18'28.50"W, 33d37'35.61"N)
Center      (  456080.000, 3735960.000) (117d28'27.39"W, 33d45'52.46"N)
Band 1 Block=512x16 Type=Byte, ColorInterp=Gray

This output may be all you need. If you want to do this programmaticly in python however, this is how you get the same info.

If the coordinate system is a PROJCS like the example above you are dealing with a projected coordinate system. A projected coordiante system is a representation of the spheroidal earth's surface, but flattened and distorted onto a plane. If you want the latitude and longitude, you need to convert the coordinates to the geographic coordinate system that you want.

Sadly, not all latitude/longitude pairs are created equal, being based upon different spheroidal models of the earth. In this example, I am converting to WGS84, the geographic coordinate system favoured in GPSs and used by all the popular web mapping sites. The coordinate system is defined by a well defined string. A catalogue of them is available from spatial ref, see for example WGS84.

from osgeo import osr, gdal

# get the existing coordinate system
ds = gdal.Open('path/to/file')
old_cs= osr.SpatialReference()
old_cs.ImportFromWkt(ds.GetProjectionRef())

# create the new coordinate system
wgs84_wkt = """
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0,
        AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.01745329251994328,
        AUTHORITY["EPSG","9122"]],
    AUTHORITY["EPSG","4326"]]"""
new_cs = osr.SpatialReference()
new_cs .ImportFromWkt(wgs84_wkt)

# create a transform object to convert between coordinate systems
transform = osr.CoordinateTransformation(old_cs,new_cs) 

#get the point to transform, pixel (0,0) in this case
width = ds.RasterXSize
height = ds.RasterYSize
gt = ds.GetGeoTransform()
minx = gt[0]
miny = gt[3] + width*gt[4] + height*gt[5] 

#get the coordinates in lat long
latlong = transform.TransformPoint(minx,miny) 

Hopefully this will do what you want.