Wie erhalte ich Rastereckkoordinaten mit Python GDAL-Bindungen?

30

Gibt es eine Möglichkeit, die Eckkoordinaten (in Grad Lat / Long) aus einer Rasterdatei mit den Python-Bindungen von gdal zu ermitteln?

Einige Online-Suchanfragen haben mich davon überzeugt, dass dies nicht der Fall ist. Daher habe ich eine Problemlösung entwickelt, indem ich die gdalinfo-Ausgabe analysiert habe. Sie ist etwas grundlegend, aber ich dachte, dass dies möglicherweise Zeit für Leute spart, die mit Python weniger vertraut sind. Es funktioniert auch nur, wenn gdalinfo die geografischen Koordinaten zusammen mit den Eckkoordinaten enthält, was meines Erachtens nicht immer der Fall ist.

Ist hier meine Problemumgehung, hat jemand bessere Lösungen?

Mit gdalinfo auf einem geeigneten Raster erhalten Sie in der Mitte der Ausgabe etwa Folgendes:

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)

Dieser Code funktioniert für Dateien, deren gdalinfo so aussieht. Ich glaube, dass die Koordinaten manchmal in Grad und Dezimalstellen anstatt in Grad, Minuten und Sekunden angegeben werden. Es sollte trivial sein, den Code für diese Situation anzupassen.

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

Und das gibt mir:

[-4.33429444 -5.82895833 -4.33429444 -5.82895833 -5.081625  ] 
[ 137.12275833  137.12203333  137.74633889  137.74706111  137.43454722]
EddyTheB
quelle
Vielleicht im Zusammenhang: gis.stackexchange.com/questions/33330/…
AndreJ

Antworten:

29

Hier ist eine andere Möglichkeit, ohne ein externes Programm aufzurufen.

Dabei werden die Koordinaten der vier Ecken von der Geotransformation abgerufen und mit osr.CoordinateTransformation in lon / lat neu projiziert.

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)

Einige Code aus dem Metageta- Projekt, osr.CoordinateTransformation Idee aus dieser Antwort

user2856
quelle
Super, danke. Und es funktioniert unabhängig davon, ob die entsprechenden Zeilen in der Ausgabe von gdalinfo vorhanden sind.
EddyTheB
Ich denke, das wird besser mit tgt_srs = src_srs.CloneGeogCS (). Meine anfänglichen Raster sind Bilder des Mars, daher ist die Verwendung von EPSG (4326) nicht ideal, aber CloneGeogCS () scheint sich nur von projizierten zu geografischen Koordinaten zu ändern.
EddyTheB
Sicher. Ich habe die Antwort aktualisiert, um CloneGeogCS zu verwenden. Ich habe jedoch nur versucht, die Verwendung der Funktionen GetExtent und ReprojectCoords zu demonstrieren. Sie können alles verwenden, was Sie als Ziel srs wollen.
user2856
Ja, danke, die hätte ich sonst nie gefunden.
EddyTheB
Was ist, wenn Sie einen Datensatz haben, der keine Projektion hat und RPCs spezifiziert? Der Import von der wkt-Funktion schlägt fehl. Es ist möglich, die Ausdehnung mit einem Transformator zu berechnen, aber ich habe mich gefragt, ob es einen Weg mit der obigen Methode gibt.
Thomas
41

Dies kann in weitaus weniger Codezeilen erfolgen

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, ulyIst die obere linke Ecke lrx, lryist die rechte untere Ecke

Die osr-Bibliothek (Teil von gdal) kann verwendet werden, um die Punkte in ein beliebiges Koordinatensystem zu transformieren. Für einen einzelnen Punkt:

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)

Um eine ganze Rasterbild zu projizieren eine weitaus komplizierte Materie wäre, aber GDAL> = 2.0 bietet eine einfache Lösung für diese auch: gdal.Warp.

James
quelle
Dies ist die pythonische Antwort für die Ausdehnung - eine ebenso pythonische Lösung für die Neuprojektion wäre fantastisch gewesen. Das heißt, ich verwende die Ergebnisse in PostGIS, übergebe also nur die nicht transformierte Ausdehnung und ST_Transform(ST_SetSRID(ST_MakeBox2D([die Ergebnisse] ),28355),4283). (Ein Streitpunkt - das 'T' in src.GetGeoTransform()sollte groß geschrieben werden).
GT.
@GT. Ein Beispiel wurde hinzugefügt
James
1

Ich habe es so gemacht ... es ist ein wenig hartcodiert, aber wenn sich nichts an der gdalinfo ändert, funktioniert es für UTM-projizierte Bilder!

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
quelle
2
Dies ist ziemlich fragil, da es darauf ankommt, dass sowohl gdalinfodie Verfügbarkeit auf dem Pfad eines Benutzers (was insbesondere unter Windows nicht immer der Fall ist) als auch eine gedruckte Ausgabe ohne strenge Schnittstelle analysiert werden, dh dass für den korrekten Abstand diese "magischen Zahlen" verwendet werden. Es ist unnötig, wenn gdal umfassende Python-Bindungen bereitstellt, die dies auf eine explizitere und robustere Weise tun können
James
1

Wenn Ihr Raster gedreht wird, wird die Berechnung etwas komplizierter, da Sie jeden der sechs affinen Transformationskoeffizienten berücksichtigen müssen. Erwägen Sie die Verwendung von affin, um eine gedrehte affine Transformation / Geotransformation zu transformieren:

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|

Jetzt können Sie die vier Eckkoordinatenpaare erzeugen:

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

Und wenn Sie auch die gitterbasierten Grenzen benötigen (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
quelle
0

Ich glaube an neuere Versionen des OSGEO / GDAL-Moduls für Python, mit dem man GDAL-Dienstprogramme direkt aus dem Code aufrufen kann, ohne dass Systemaufrufe erforderlich sind. Beispiel: Anstatt einen Unterprozess zum Aufrufen zu verwenden:

gdalinfo kann mit gdal.Info (der_Dateiname) aufgerufen werden, um die Metadaten / Anmerkungen der Datei anzuzeigen

Anstatt den Unterprozess zum Aufrufen von: gdalwarp zu verwenden, können Sie gdal.Warp aufrufen, um das Warping für ein Raster auszuführen.

Die Liste der derzeit als interne Funktion verfügbaren GDAL-Dienstprogramme: http://erouault.blogspot.com/2015/10/gdal-and-ogr-utilities-as-library.html

Sinooshka
quelle