Obtenir la latitude et la longitude à partir d'un fichier GeoTIFF

Question

En utilisant GDAL en Python, comment obtenez-vous la latitude et la longitude d'un fichier GeoTIFF?

Les GeoTIFF ne semblent pas stocker d'informations de coordonnées. Au lieu de cela, ils stockent les coordonnées XY Origin. Cependant, les coordonnées XY ne fournissent pas la latitude et la longitude du coin supérieur gauche et du coin inférieur gauche.

Il semble que j'aurai besoin de faire quelques calculs pour résoudre ce problème, mais je n'ai pas la moindre idée de l'endroit où commencer.

Quelle procédure est nécessaire pour effectuer cette opération?

Je sais que la méthode GetGeoTransform() est importante pour cela, mais je ne sais pas quoi faire avec.

Answer 1

Pour obtenir les coordonnées des coins de votre geotiff, procédez comme suit:

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] 

Cependant, ceux-ci pourraient ne pas être au format latitude/longitude. Comme Justin l'a noté, votre géotiff sera stocké avec une sorte de système de coordonnées. Si vous ne savez pas ce que le système de coordonnées est, vous pouvez trouver en exécutant gdalinfo:

gdalinfo ~/somedir/somefile.tif 

qui sort:

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 

Cette sortie peut être tout ce dont vous avez besoin. Si vous voulez le faire par programmation en python, voici comment vous obtenez la même information.

Si le système de coordonnées est un PROJCS comme dans l'exemple ci-dessus, vous avez affaire à un système de coordonnées projeté. Un système de coordonnées projeté est un representation of the spheroidal earth's surface, but flattened and distorted onto a plane. Si vous voulez la latitude et la longitude, vous devez convertir les coordonnées au système de coordonnées géographiques que vous voulez. Malheureusement, toutes les paires latitude/longitude ne sont pas égales, étant basées sur différents modèles sphéroïdaux de la Terre. Dans cet exemple, je convertis en WGS84, le système de coordonnées géographiques privilégié dans les GPS et utilisé par tous les sites de cartographie Web populaires. Le système de coordonnées est défini par une chaîne bien définie. Un catalogue d'entre eux est disponible à partir de spatial ref, voir par exemple 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) 

Espérons que cela fera ce que vous voulez.

Answer 2

Je ne sais pas si cela est une réponse complète, mais this site dit:

Le x/y dimensions de carte sont appelés et Ordonnée Abscisse. Pour les ensembles de données dans un système de coordonnées géographiques, ceux-ci pourraient contenir la longitude et la latitude. Pour les systèmes de coordonnées projetés, ils sont normalement l'est et le nord dans le système de coordonnées projeté. Pour les images non géoréférencées, l'est et le nord seraient simplement les décalages pixel/ligne de chaque pixel (comme cela est implicite par une géotransformation d'unité).

donc ils peuvent être réellement la longitude et la latitude.