¿Cómo obtener coordenadas de esquina ráster usando enlaces Python GDAL?

30

¿Hay alguna manera de obtener las coordenadas de las esquinas (en grados lat / long) de un archivo ráster usando los enlaces de Python de gdal?

Algunas búsquedas en línea me han convencido de que no existe, así que he desarrollado una solución analizando el resultado de gdalinfo, es algo básico, pero pensé que podría ahorrar algo de tiempo para las personas que podrían estar menos cómodas con Python. También solo funciona si gdalinfo contiene las coordenadas geográficas junto con las coordenadas de las esquinas, lo cual no creo que sea siempre el caso.

Aquí está mi solución, ¿alguien tiene alguna solución mejor?

gdalinfo en un ráster apropiado da como resultado algo así a la mitad de la salida:

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)

Este código funcionará en archivos cuyo gdalinfo se vea así. Creo que a veces las coordenadas estarán en grados y decimales, en lugar de grados, minutos y segundos; Debería ser trivial ajustar el código para esa situación.

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

Y eso me da:

[-4.33429444 -5.82895833 -4.33429444 -5.82895833 -5.081625  ] 
[ 137.12275833  137.12203333  137.74633889  137.74706111  137.43454722]
EddyTheB
fuente
Quizás relacionado: gis.stackexchange.com/questions/33330/…
AndreJ

Respuestas:

29

Aquí hay otra forma de hacerlo sin llamar a un programa externo.

Lo que esto hace es obtener las coordenadas de las cuatro esquinas de la geotransformación y reproyectarlas a lon / lat usando 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)

Parte del código de la metageta proyecto, idea osr.CoordinateTransformation de esta respuesta

usuario2856
fuente
Genial gracias. Y funciona independientemente de si existen las líneas apropiadas en la salida de gdalinfo.
EddyTheB
Creo que esto será mejor con tgt_srs = src_srs.CloneGeogCS (). Mis rásteres iniciales son imágenes de Marte, por lo que usar EPSG (4326) no es ideal, pero CloneGeogCS () parece cambiar de coordenadas proyectadas a geográficas.
EddyTheB
Sin lugar a duda. He actualizado la respuesta para usar CloneGeogCS. Sin embargo, solo estaba tratando de demostrar el uso de las funciones GetExtent y ReprojectCoords. Puedes usar lo que quieras como objetivo srs.
usuario2856
Sí, gracias, nunca los habría encontrado de otra manera.
EddyTheB
¿Qué sucede si tiene un conjunto de datos que no tiene proyección y especifica RPC? La importación desde la función wkt falla. Es posible calcular la extensión usando un transformador, pero me preguntaba si había alguna manera con el método anterior.
Thomas
41

Esto se puede hacer en muchas menos líneas de código

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, ulyes la esquina superior izquierda lrx, lryes la esquina inferior derecha

La biblioteca osr (parte de gdal) se puede usar para transformar los puntos en cualquier sistema de coordenadas. Por un solo punto:

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)

Reproyectar una imagen de trama entera sería un asunto mucho más complicado, pero GDAL> = 2.0 ofrece una solución fácil para esto también: gdal.Warp.

James
fuente
Esta es la respuesta Pythonic para la extensión: una solución igualmente Pythonic para la reproyección habría sido increíble. Dicho esto, utilizo los resultados en PostGIS, así que solo paso la extensión no transformada y ST_Transform(ST_SetSRID(ST_MakeBox2D([los resultados] ),28355),4283). (Una objeción: la 'T' src.GetGeoTransform()debe estar en mayúscula).
GT.
@GT. Se agregó un ejemplo
James
1

Lo he hecho de esta manera ... está un poco codificado, pero si nada cambia con el gdalinfo, ¡funcionará para las imágenes proyectadas de UTM!

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]
Emiliano
fuente
2
Esto es bastante frágil, ya que depende tanto de gdalinfoestar disponible en la ruta del usuario (no siempre es el caso, especialmente en Windows) como de analizar una salida impresa que no tiene una interfaz estricta, es decir, confiar en esos 'números mágicos' para un espaciado correcto. No es necesario cuando gdal proporciona enlaces completos de python que pueden hacer esto de una manera más explícita y robusta
James
1

Si su ráster se gira, entonces las matemáticas se vuelven un poco más complicadas, ya que debe considerar cada uno de los seis coeficientes de transformación afines. Considere usar afín para transformar una transformación afín rotada / geotransformada:

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|

Ahora puede generar los cuatro pares de coordenadas de esquina:

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

Y si también necesita los límites basados ​​en la cuadrícula (xmin, ymin, xmax, ymax):

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

Creo que las versiones más recientes del módulo OSGEO / GDAL para python one pueden llamar directamente a las utilidades GDAL desde el código sin involucrar llamadas al sistema. por ejemplo, en lugar de usar un subproceso para llamar:

gdalinfo se puede llamar a gdal.Info (the_name_of_the_file) para tener una exposición de los metadatos / anotaciones del archivo

o en lugar de usar un subproceso para llamar: gdalwarp se puede llamar a gdal. Warp para realizar la deformación en una trama.

La lista de utilidades GDAL actualmente disponibles como una función interna: http://erouault.blogspot.com/2015/10/gdal-and-ogr-utilities-as-library.html

Sinooshka
fuente