TransWikia.com

How to get raster corner coordinates using Python GDAL bindings?

Geographic Information Systems Asked by EddyTheB on March 30, 2021

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]

6 Answers

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

Correct answer by user2856 on March 30, 2021

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]

Answered by Emiliano on March 30, 2021

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.

Answered by James on March 30, 2021

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)

Answered by Mike T on March 30, 2021

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

Answered by Sinooshka on March 30, 2021

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

Answered by Francisco Javier López Andreu on March 30, 2021

Add your own answers!

Ask a Question

Get help from others!

© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP