How to get raster corner coordinates using Python GDAL bindings?
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
.
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(gt,cols,rows):
''' Return list of corner coordinates from a geotransform
@type gt: C{tuple/list}
@param gt: geotransform
@type cols: C{int}
@param cols: number of columns in the dataset
@type rows: C{int}
@param rows: number of rows in the dataset
@rtype: C{[float,...,float]}
@return: coordinates of each corner
'''
ext=[]
xarr=[0,cols]
yarr=[0,rows]
for px in xarr:
for py in yarr:
x=gt[0]+(px*gt[1])+(py*gt[2])
y=gt[3]+(px*gt[4])+(py*gt[5])
ext.append([x,y])
print x,y
yarr.reverse()
return ext
def ReprojectCoords(coords,src_srs,tgt_srs):
''' Reproject a list of x,y coordinates.
@type geom: C{tuple/list}
@param geom: List of [[x,y],...[x,y]] coordinates
@type src_srs: C{osr.SpatialReference}
@param src_srs: OSR SpatialReference object
@type tgt_srs: C{osr.SpatialReference}
@param tgt_srs: OSR SpatialReference object
@rtype: C{tuple/list}
@return: List of transformed [[x,y],...[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)
gt=ds.GetGeoTransform()
cols = ds.RasterXSize
rows = ds.RasterYSize
ext=GetExtent(gt,cols,rows)
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 the metageta project, osr.CoordinateTransformation idea from this answer
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