43

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: break 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]
Improve this question
1

7 Answers 7

Reset to default
73

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() source.ImportFromWkt(src.GetProjection()) # The target projection target = osr.SpatialReference() target.ImportFromEPSG(4326) # 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.

Improve this answer
3
  • 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). Commented Feb 22, 2017 at 0:59
  • @GT. Added an example Commented 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 Commented Jan 26, 2022 at 21:02
41

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. """ trans_coords=[] transform = osr.CoordinateTransformation( src_srs, tgt_srs) for x,y in coords: x,y,z = transform.TransformPoint(x,y) trans_coords.append([x,y]) return trans_coords raster=r'somerasterfile.tif' ds=gdal.Open(raster) ext=GetExtent(ds) src_srs=osr.SpatialReference() src_srs.ImportFromWkt(ds.GetProjection()) #tgt_srs=osr.SpatialReference() #tgt_srs.ImportFromEPSG(4326) tgt_srs = src_srs.CloneGeogCS() geo_ext=ReprojectCoords(ext, src_srs, tgt_srs)

Some code from this answer

Improve this answer
8
  • Awesome, thanks. And it works regardless of whether the appropriate lines exist in the gdalinfo output. Commented 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. Commented 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. Commented Apr 19, 2013 at 1:17
  • Yes, thank you, I'd have never found those otherwise. Commented 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? Commented Oct 23, 2018 at 21:12
10

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') print(info['wgs84Extent']) # {'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

Improve this answer
4

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) print(transform) # | 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)
Improve this answer
1
4

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
Improve this answer
2

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]
Improve this answer
1
  • 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 Commented Feb 22, 2017 at 7:44
0

With GDAL 3.12, this can be done via any of:

ds = gdal.Open('data.tif') # extent in native projection ds.GetExtent() # extent in some other projection (WGS84 lat/lon) ds.GetExtent(osr.SpatialReference(epsg=4326)) # extent in WGS84 lon/lat ds.GetExtentWGS84LongLat()

The extent is returned as a tuple of xmin, xmax, ymin, ymax.

Improve this answer

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.