Is there a way to get the corner coordinates (in degrees lat/long) from a raster file using gdal's Python bindings?

A few searches online have convinced me that there is not, so I have developed a work around by parsing the gdalinfo output, it's somewhat basic but I thought it might save some time for people who might be less comfortable with python. It also only works if gdalinfo contains the geographic coordinates along with the corner coordinates, which I don't believe is always the case.

Here's my workaround, does anyone have any better solutions?

gdalinfo on an appropriate raster results in something like this midway through the output:

Corner Coordinates:
Upper Left  (  -18449.521, -256913.934) (137d 7'21.93"E,  4d20'3.46"S)
Lower Left  (  -18449.521, -345509.683) (137d 7'19.32"E,  5d49'44.25"S)
Upper Right (   18407.241, -256913.934) (137d44'46.82"E,  4d20'3.46"S)
Lower Right (   18407.241, -345509.683) (137d44'49.42"E,  5d49'44.25"S)
Center      (     -21.140, -301211.809) (137d26'4.37"E,  5d 4'53.85"S)

This code will work on files who's gdalinfo look like that. I believe sometimes the coordinates will be in degrees and decimals, rather than degrees, minutes and seconds; it ought to be trivial to adjust the code for that situation.

import numpy as np
import subprocess

def GetCornerCoordinates(FileName):
    GdalInfo = subprocess.check_output('gdalinfo {}'.format(FileName), shell=True)
    GdalInfo = GdalInfo.split('/n') # Creates a line by line list.
    CornerLats, CornerLons = np.zeros(5), np.zeros(5)
    GotUL, GotUR, GotLL, GotLR, GotC = False, False, False, False, False
    for line in GdalInfo:
        if line[:10] == 'Upper Left':
            CornerLats[0], CornerLons[0] = GetLatLon(line)
            GotUL = True
        if line[:10] == 'Lower Left':
            CornerLats[1], CornerLons[1] = GetLatLon(line)
            GotLL = True
        if line[:11] == 'Upper Right':
            CornerLats[2], CornerLons[2] = GetLatLon(line)
            GotUR = True
        if line[:11] == 'Lower Right':
            CornerLats[3], CornerLons[3] = GetLatLon(line)
            GotLR = True
        if line[:6] == 'Center':
            CornerLats[4], CornerLons[4] = GetLatLon(line)
            GotC = True
        if GotUL and GotUR and GotLL and GotLR and GotC:
    return CornerLats, CornerLons 

def GetLatLon(line):
    coords = line.split(') (')[1]
    coords = coords[:-1]
    LonStr, LatStr = coords.split(',')
    # Longitude
    LonStr = LonStr.split('d')    # Get the degrees, and the rest
    LonD = int(LonStr[0])
    LonStr = LonStr[1].split('\'')# Get the arc-m, and the rest
    LonM = int(LonStr[0])
    LonStr = LonStr[1].split('"') # Get the arc-s, and the rest
    LonS = float(LonStr[0])
    Lon = LonD + LonM/60. + LonS/3600.
    if LonStr[1] in ['W', 'w']:
        Lon = -1*Lon
    # Same for Latitude
    LatStr = LatStr.split('d')
    LatD = int(LatStr[0])
    LatStr = LatStr[1].split('\'')
    LatM = int(LatStr[0])
    LatStr = LatStr[1].split('"')
    LatS = float(LatStr[0])
    Lat = LatD + LatM/60. + LatS/3600.
    if LatStr[1] in ['S', 's']:
        Lat = -1*Lat
    return Lat, Lon

FileName = Image.cub
# Mine's an ISIS3 cube file.
CornerLats, CornerLons = GetCornerCoordinates(FileName)
# UpperLeft, LowerLeft, UpperRight, LowerRight, Centre
print CornerLats
print CornerLons

And that gives me:

[-4.33429444 -5.82895833 -4.33429444 -5.82895833 -5.081625  ] 
[ 137.12275833  137.12203333  137.74633889  137.74706111  137.43454722]

6 Answers 6


This can be done in far fewer lines of code

src = gdal.Open(path goes here)
ulx, xres, xskew, uly, yskew, yres  = src.GetGeoTransform()
lrx = ulx + (src.RasterXSize * xres)
lry = uly + (src.RasterYSize * yres)

ulx, uly is the upper left corner, lrx, lry is the lower right corner

The osr library (part of gdal) can be used to transform the points to any coordinate system. For a single point:

from osgeo import ogr
from osgeo import osr

# Setup the source projection - you can also import from epsg, proj4...
source = osr.SpatialReference()

# The target projection
target = osr.SpatialReference()

# Create the transform - this can be used repeatedly
transform = osr.CoordinateTransformation(source, target)

# Transform the point. You can also create an ogr geometry and use the more generic `point.Transform()`
transform.TransformPoint(ulx, uly)

To reproject a whole raster image would be a far more complicated matter, but GDAL >= 2.0 offers an easy solution for this too: gdal.Warp.

  • This is the Pythonic answer for the extent - an equally-Pythonic solution for the reprojection would have been awesome, That said - I use the results in PostGIS, so I just pass the un-transformed extent and ST_Transform(ST_SetSRID(ST_MakeBox2D( [the results] ),28355),4283) it. (One quibble - the 'T' in src.GetGeoTransform() should be capitalised).
    – GT.
    Feb 22, 2017 at 0:59
  • @GT. Added an example
    – James
    Feb 22, 2017 at 7:37
  • you need to do ul_lat, ul_lon, _ = transform.TransformPoint(ulx, uly), as python doesn't modify ulx and uly in place Jan 26, 2022 at 21:02

Here's another way to do it without calling an external program.

What this does is get the coordinates of the four corners from the geotransform and reproject them to lon/lat using osr.CoordinateTransformation.

from osgeo import gdal,ogr,osr

def GetExtent(ds):
    """ Return list of corner coordinates from a gdal Dataset """
    xmin, xpixel, _, ymax, _, ypixel = ds.GetGeoTransform()
    width, height = ds.RasterXSize, ds.RasterYSize
    xmax = xmin + width * xpixel
    ymin = ymax + height * ypixel

    return (xmin, ymax), (xmax, ymax), (xmax, ymin), (xmin, ymin)

def ReprojectCoords(coords,src_srs,tgt_srs):
    """ Reproject a list of x,y coordinates. """
    transform = osr.CoordinateTransformation( src_srs, tgt_srs)
    for x,y in coords:
        x,y,z = transform.TransformPoint(x,y)
    return trans_coords



tgt_srs = src_srs.CloneGeogCS()

geo_ext=ReprojectCoords(ext, src_srs, tgt_srs)

Some code from this answer

  • Awesome, thanks. And it works regardless of whether the appropriate lines exist in the gdalinfo output.
    – EddyTheB
    Apr 16, 2013 at 18:01
  • I think this will be better with tgt_srs = src_srs.CloneGeogCS(). My initial rasters are images of Mars, so using EPSG(4326) isn't ideal, but CloneGeogCS() appears to just change from projected to geographic coordinates.
    – EddyTheB
    Apr 19, 2013 at 0:49
  • For sure. I've updated the answer to use CloneGeogCS. However, I was just trying to demonstrate the use of the GetExtent and ReprojectCoords functions. You can use anything you want as target srs.
    – user2856
    Apr 19, 2013 at 1:17
  • Yes, thank you, I'd have never found those otherwise.
    – EddyTheB
    Apr 19, 2013 at 2:28
  • What if you have a dataset which has no projection and specifies RPCs? The import from wkt function fails. It is possible to calculate the extent using a transformer but I was wondering if there was a way with the method above?
    – Thomas
    Oct 23, 2018 at 21:12

I believe in more recent versions of OSGEO/GDAL module for python one can directly call GDAL utilities from code without involving system calls. for example instead of using subprocess to call :

gdalinfo one can call gdal.Info(the_name_of_the_file) to have an exposure of the file metadata/annotation, like this:

    from osgeo import gdal
    raster_path = '....tif'
    info = gdal.Info(raster_path, format='json')
    # {'type': 'Polygon',
    #  'coordinates': [[[6.0, 44.0], [6.0, 44.0], [7.0, 4.0], [7.0, 44.0], [6.0, 44.0]]]}

or instead of using subprocess to call: gdalwarp one can cal gdal.Warp to do the warping on a raster.

The list of GDAL utilities currently available as an internal function : http://erouault.blogspot.com/2015/10/gdal-and-ogr-utilities-as-library.html


If your raster is rotated, then the math gets a bit more complicated, as you need to consider each of the six affine transformation coefficients. Consider using affine to transform a rotated affine transformation / geotransform:

from affine import Affine

# E.g., from a GDAL DataSet object:
# gt = ds.GetGeoTransform()
# ncol = ds.RasterXSize
# nrow = ds.RasterYSize

# or to work with a minimal example
gt = (100.0, 17.320508075688775, 5.0, 200.0, 10.0, -8.660254037844387)
ncol = 10
nrow = 15

transform = Affine.from_gdal(*gt)
# | 17.32, 5.00, 100.00|
# | 10.00,-8.66, 200.00|
# | 0.00, 0.00, 1.00|

Now you can generate the four corner coordinate pairs:

c0x, c0y = transform.c, transform.f  # upper left
c1x, c1y = transform * (0, nrow)     # lower left
c2x, c2y = transform * (ncol, nrow)  # lower right
c3x, c3y = transform * (ncol, 0)     # upper right

And if you also need the grid-based bounds (xmin, ymin, xmax, ymax):

xs = (c0x, c1x, c2x, c3x)
ys = (c0y, c1y, c2y, c3y)
bounds = min(xs), min(ys), max(xs), max(ys)

With rasterio is very simple, I've done it like this:

ds_raster = rasterio.open(raster_path)
bounds = ds_raster.bounds
left= bounds.left
bottom = bounds.bottom
right = bounds.right
top = bounds.top

I've done this way... it is a little hard-coded but if nothing changes with the gdalinfo, it will work for UTM projected images!

imagefile= /pathto/image
p= subprocess.Popen(["gdalinfo", "%s"%imagefile], stdout=subprocess.PIPE)
out,err= p.communicate()
ul= out[out.find("Upper Left")+15:out.find("Upper Left")+38]
lr= out[out.find("Lower Right")+15:out.find("Lower Right")+38]
  • 2
    This is fairly fragile as it relies on both gdalinfo being available on a users' path (not always the case, especially on windows) and parsing a printed output which has no strict interface - i.e. relying on those 'magic numbers' for correct spacing. It is unnecessary when gdal provides comprehensive python bindings which can do this in a more explicit and robust manner
    – James
    Feb 22, 2017 at 7:44

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.