Latitude und Longitude aus einer GeoTIFF-Datei beziehen

Question

Wie erhalten Sie mit GDAL in Python den Längen- und Breitengrad einer GeoTIFF-Datei?

GeoTIFFs scheinen keine Koordinateninformationen zu speichern. Stattdessen speichern sie die XY-Koordinaten. Die XY-Koordinaten liefern jedoch nicht den Breiten- und Längengrad der oberen linken Ecke und der linken unteren Ecke.

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

Welche Vorgehensweise ist erforderlich, damit dies durchgeführt wird?

Ich weiß, dass die GetGeoTransform() Methode dafür wichtig ist, aber ich weiß nicht, was damit von dort zu tun ist.

Answer 1

Um die Koordinaten der Ecken Ihres geotiff gehen Sie wie folgt zu erhalten:

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] 

jedoch diese möglicherweise nicht in Breite/Länge-Format vorliegen. Wie Justin bemerkte, wird dein GeoTiff mit einer Art Koordinatensystem gespeichert. Wenn Sie nicht wissen, welches Koordinatensystem ist, können Sie herausfinden, gdalinfo durch Ausführen von:

gdalinfo ~/somedir/somefile.tif 

Welche Ausgänge:

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 

Dieser Ausgang alles was Sie brauchen werden können. Wenn Sie dies programmatisch in Python tun möchten, erhalten Sie auf diese Weise die gleichen Informationen.

Wenn das Koordinatensystem wie im obigen Beispiel ein PROJCS ist, handelt es sich um ein projiziertes Koordinatensystem. Ein projiziertes Koordinatensystem ist ein representation of the spheroidal earth's surface, but flattened and distorted onto a plane. Wenn Sie den Längen- und Breitengrad möchten, müssen Sie die Koordinaten in das gewünschte geografische Koordinatensystem konvertieren.

Leider sind nicht alle Breiten-/Längengradpaare gleich, da sie auf verschiedenen Kugelmodellen der Erde basieren. In diesem Beispiel konvertiere ich nach WGS84, dem geografischen Koordinatensystem, das in GPS-Geräten bevorzugt wird und von allen gängigen Web-Mapping-Sites verwendet wird. Das Koordinatensystem wird durch eine gut definierte Zeichenkette definiert. Ein Katalog davon ist erhältlich von spatial ref, 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(x,y) 

Hoffentlich wird dies tun, was Sie wollen.

Answer 2

Ich weiß nicht, ob dies eine vollständige Antwort, aber this site sagt:

Die x/y Karte Dimensionen Rechts- und Hochwert bezeichnet werden. Für Datensätze in einem geographischen Koordinatensystem würden diese den Längen- und Breitengrad enthalten. Für projizierte Koordinatensysteme wären sie normalerweise der Ost- und Nordwert im projizierten Koordinatensystem. Bei nicht georeferenzierten Bildern wären der Ost- und der Nordwert nur die Pixel-/Linien-Offsets jedes Pixels (was durch eine Einheits-Geotransformation impliziert wird).

so dass sie tatsächlich Länge und Breite sein können.