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]
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
Get help from others!
Recent Questions
Recent Answers
© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP