Obtener latitud y longitud de un archivo GeoTIFF

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

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

Parece que tendré que hacer algunos cálculos matemáticos para resolver este problema, pero no tengo 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í.

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, estos pueden no estar en formato de latitud / longitud. Como señaló Justin, su geotiff se almacenará con algún tipo de sistema de coordenadas. Si no sabe qué sistema de coordenadas es, puede averiguarlo 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 

Esta salida puede ser todo lo que necesita. Sin embargo, si quieres hacer esto programáticamente en Python, esta es la forma en que obtienes la misma información.

Si el sistema de coordenadas es un PROJCS como en el ejemplo anterior, se trata de un sistema de coordenadas proyectado. Un sistema de coordenadas proyectado es una representación de la superficie de la tierra esferoidal, pero aplanada y distorsionada en un plano . 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 son creados iguales, ya que se basan 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 de mapas web populares. El sistema de coordenadas está definido por una cadena bien definida. Un catálogo de ellos está disponible en la 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) 

Esperemos que esto haga lo que quieras.

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

Las dimensiones del mapa x / y se denominan este y norte. Para los conjuntos de datos en un sistema de coordenadas geográficas, estos tendrían la longitud y la latitud. Para los sistemas de coordenadas proyectadas, normalmente serían el este y el norte en el sistema de coordenadas proyectadas. Para imágenes no referenciadas, el este y el norte serían solo las compensaciones de píxel / línea de cada píxel (como lo implica una geotransforma de unidad).

así que en realidad pueden ser longitud y latitud.