Abrufen von Breiten- und Längengraden aus einer GeoTIFF-Datei

Question

Wie kann man mit GDAL in Python den Breiten- und Längengrad einer GeoTIFF-Datei ermitteln?

GeoTIFFs scheinen keine Koordinateninformationen zu speichern. Stattdessen speichern sie die XY-Ursprungskoordinaten. Die XY-Koordinaten geben jedoch nicht den Längen- und Breitengrad der linken oberen Ecke und der linken unteren Ecke an.

Es scheint, dass ich etwas rechnen muss, um dieses Problem zu lösen, aber ich habe keine Ahnung, wo ich anfangen soll.

Welches Verfahren ist für die Durchführung dieser Maßnahme erforderlich?

Ich weiß, dass die GetGeoTransform() Methode ist dafür wichtig, aber ich weiß nicht, was ich damit anfangen soll.

Answer 1

Um die Koordinaten der Ecken Ihres Geotiffs zu erhalten, gehen Sie wie folgt vor:

from osgeo import gdal
ds = gdal.Open('path/to/file')
width = ds.RasterXSize
height = ds.RasterYSize
gt = ds.GetGeoTransform()
minx = gt[0]
miny = gt[3] + width*gt[4] + height*gt[5] 
maxx = gt[0] + width*gt[1] + height*gt[2]
maxy = gt[3] 

Diese sind jedoch möglicherweise nicht im Format Breitengrad/Längengrad. Wie Justin bemerkte, wird Ihr Geotiff mit einer Art Koordinatensystem gespeichert. Wenn Sie nicht wissen, um welches Koordinatensystem es sich handelt, können Sie es herausfinden, indem Sie gdalinfo :

gdalinfo ~/somedir/somefile.tif 

Welche Ausgaben:

Driver: GTiff/GeoTIFF
Size is 512, 512
Coordinate System is:
PROJCS["NAD27 / UTM zone 11N",
    GEOGCS["NAD27",
        DATUM["North_American_Datum_1927",
            SPHEROID["Clarke 1866",6378206.4,294.978698213901]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",-117],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1]]
Origin = (440720.000000,3751320.000000)
Pixel Size = (60.000000,-60.000000)
Corner Coordinates:
Upper Left  (  440720.000, 3751320.000) (117d38'28.21"W, 33d54'8.47"N)
Lower Left  (  440720.000, 3720600.000) (117d38'20.79"W, 33d37'31.04"N)
Upper Right (  471440.000, 3751320.000) (117d18'32.07"W, 33d54'13.08"N)
Lower Right (  471440.000, 3720600.000) (117d18'28.50"W, 33d37'35.61"N)
Center      (  456080.000, 3735960.000) (117d28'27.39"W, 33d45'52.46"N)
Band 1 Block=512x16 Type=Byte, ColorInterp=Gray

Diese Ausgabe könnte alles sein, was Sie brauchen. Wenn Sie dies jedoch programmatisch in Python tun möchten, erhalten Sie die gleichen Informationen auf diese Weise.

Wenn das Koordinatensystem ein PROJCS wie im obigen Beispiel haben Sie es mit einem projizierten Koordinatensystem zu tun. Ein projiziertes Koordiantensystem ist ein Darstellung der kugelförmigen Erdoberfläche, jedoch abgeflacht und verzerrt auf einer Ebene . Wenn Sie den Breiten- und Längengrad benötigen, müssen Sie die Koordinaten in das gewünschte geografische Koordinatensystem umwandeln.

Leider sind nicht alle Breiten-/Längenpaare gleich, da sie auf unterschiedlichen Sphäroidmodellen der Erde beruhen. In diesem Beispiel konvertiere ich in WGS84 das geografische Koordinatensystem, das in GPS-Geräten bevorzugt und von allen beliebten Web-Mapping-Seiten verwendet wird. Das Koordinatensystem wird durch eine genau definierte Zeichenfolge festgelegt. Ein entsprechender Katalog ist erhältlich bei räumlicher Bezug siehe zum Beispiel WGS84 .

from osgeo import osr, gdal

# get the existing coordinate system
ds = gdal.Open('path/to/file')
old_cs= osr.SpatialReference()
old_cs.ImportFromWkt(ds.GetProjectionRef())

# create the new coordinate system
wgs84_wkt = """
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0,
        AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.01745329251994328,
        AUTHORITY["EPSG","9122"]],
    AUTHORITY["EPSG","4326"]]"""
new_cs = osr.SpatialReference()
new_cs .ImportFromWkt(wgs84_wkt)

# create a transform object to convert between coordinate systems
transform = osr.CoordinateTransformation(old_cs,new_cs) 

#get the point to transform, pixel (0,0) in this case
width = ds.RasterXSize
height = ds.RasterYSize
gt = ds.GetGeoTransform()
minx = gt[0]
miny = gt[3] + width*gt[4] + height*gt[5] 

#get the coordinates in lat long
latlong = transform.TransformPoint(minx,miny) 

Ich hoffe, dass dies Ihren Wünschen entspricht.

Answer 2

Ich weiß nicht, ob dies eine vollständige Antwort ist, aber diese Seite dice:

Die x/y-Dimensionen der Karte werden als Ost- und Nordrichtung bezeichnet. Bei Datensätzen in einem geografischen Koordinatensystem würden diese den Längen- und Breitengrad enthalten. Bei projizierten Koordinatensystemen handelt es sich in der Regel um den Ost- und Nordwert im projizierten Koordinatensystem. Bei nicht referenzierten Bildern wären Ost- und Nordrichtung einfach die Pixel-/Linienabstände der einzelnen Pixel (wie bei einer einheitlichen Geotransformation).

Es kann sich also tatsächlich um Längen- und Breitengrade handeln.