Obtener latitud y longitud desde un archivo GeoTIFF

Question

Usando GDAL en Python, ¿cómo se obtiene la latitud y la longitud de un archivo GeoTIFF?

Los GeoTIFF no parecen almacenar ninguna información de coordenadas. En cambio, almacenan las coordenadas XY Origen. Sin embargo, las coordenadas XY no proporcionan la latitud y la longitud de la esquina superior izquierda y la esquina inferior izquierda.

Parece que tendré que hacer algunos cálculos para resolver este problema, pero no tengo ni idea de por dónde empezar.

¿Qué procedimiento se requiere para realizar esto?

Sé que el método GetGeoTransform() es importante para esto, sin embargo, no sé qué hacer con él desde allí.

Answer 1

Para obtener las coordenadas de las esquinas de su geotiff haga lo siguiente:

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] 

Sin embargo, éstos no podrían estar en formato de latitud/longitud. Como notó Justin, tu geotiff se almacenará con algún tipo de sistema de coordenadas. Si usted no sabe qué sistema de coordenadas es, se puede averiguar mediante la ejecución gdalinfo:

gdalinfo ~/somedir/somefile.tif 

que da salida:

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 

Esta salida puede ser todo lo que necesita. Sin embargo, si quieres hacer esto programáticamente en python, así es como obtienes la misma información.

Si el sistema de coordenadas es PROJCS como en el ejemplo anterior, se trata de un sistema de coordenadas proyectadas. Un sistema de coordinación proyectado es un representation of the spheroidal earth's surface, but flattened and distorted onto a plane. Si desea la latitud y la longitud, debe convertir las coordenadas al sistema de coordenadas geográficas que desee.

Lamentablemente, no todos los pares de latitud/longitud se crean iguales, basándose en diferentes modelos esferoidales de la tierra. En este ejemplo, estoy convirtiendo a WGS84, el sistema de coordenadas geográficas preferido en los GPS y utilizado por todos los sitios populares de mapeo web. El sistema de coordenadas está definido por una cadena bien definida. Un catálogo de ellos está disponible en spatial ref, ver por ejemplo 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) 

Espero que esto haga lo que quiera.

Answer 2

No sé si esto es una respuesta completa, pero this site dice:

El x/y dimensiones mapa se llaman este y norte. Para conjuntos de datos en un sistema de coordenadas geográficas, estos mantendrían la longitud y la latitud. Para sistemas de coordenadas proyectados, normalmente serían el este y el norte en el sistema de coordenadas proyectadas. Para las imágenes no georreferenciadas, el este y el norte serían simplemente los desplazamientos de píxel/línea de cada píxel (como lo implica un geotransformado de unidad).

por lo que en realidad puede ser longitud y latitud.