python - Obtener latitud y longitud desde un archivo GeoTIFF

math geolocation

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, desde ese momento no sé qué hacer con él.

Answer 1

No sé si esta es una respuesta completa, pero este sitio dice:

Las dimensiones del mapa x / y se llaman easting y northing. 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 pueden ser de longitud y latitud.

Answer 2

Para obtener las coordenadas de las esquinas de tu geotiff haz 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, estos pueden no estar en formato de latitud / longitud. Como notó Justin, tu geotiff se almacenará con algún tipo de sistema de coordenadas. Si no sabes qué sistema de coordenadas es, puedes averiguar ejecutando gdalinfo :

gdalinfo ~/somedir/somefile.tif 

Qué salidas:

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

Este resultado 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 un PROJCS como el ejemplo anterior, se trata de un sistema de coordenadas proyectadas. Un sistema coordinado proyectado es una representación de la superficie de la tierra esferoidal, pero aplanada y distorsionada sobre un plano . Si desea la latitud y la longitud, debe convertir las coordenadas al sistema de coordenadas geográficas que desee.

Tristemente, no todos los pares de latitud / longitud son creados iguales, basándose en diferentes modelos esferoidales de la tierra. En este ejemplo, me 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 referencia espacial , 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 quieras.