Comment convertir des coordonnées affines en lat / lng?

Je suis extrêmement nouveau dans le SIG.

J'utilise gdalpour lire une carte d'utilisation des terres / couverture terrestre et j'ai besoin de choisir le lat / lng de certains types de couverture terrestre à indexer dans un ensemble de données différent qui n'est exprimé qu'en lat / lng. Malheureusement, je ne comprends pas la forme des coordonnées x et y qui me sont données par la géotransformation, en particulier la originXet originYci - dessous:

geotransform = dataset.GetGeoTransform()
originX = geotransform[0]
originY = geotransform[3]

L'impression de ces valeurs me donne des coordonnées comme (447466.693808, 4952570.40529). Quel est leur rapport avec la latitude et la longitude d'origine?

Éditer:

Voici un exemple simple de python qui m'a donné ce que je cherchais:

srs = osr.SpatialReference()
srs.ImportFromWkt(dataset.GetProjection())

srsLatLong = srs.CloneGeogCS()
ct = osr.CoordinateTransformation(srs,srsLatLong)
print ct.TransformPoint(originX,originY)

Volé à: tolatlong.py

gdal_translate va reprojeter vos données de n'importe quelle projection dans laquelle elles se trouvent vers une autre (dans ce cas, vous voulez EPSG: 4326) en utilisant:

gdal_translate -a_srs epsg:4326 srcfile new_file 

ou vous pouvez utiliser gdaltrasform pour convertir les points (et je suis sûr que vous pouvez aussi y accéder depuis Python (?))

La géotransformation est documentée sur https://gdal.org/user/raster_data_model.html . L'idée est que vous prenez les coordonnées (x, y) directement à partir du jeu de données, appliquez une transformation linéaire pour obtenir (u, v) avec

u = a*x + b*y
v = c*x + d*y

(vous pouvez prendre cela comme la définition d'une transformation linéaire), puis déplacez l'origine en ajoutant la géotransformation [0] à u et la géotransformation [3] à v. Cela donne la "transformation affine" de (x, y). Il est vraiment destiné à faire pivoter, changer d'échelle, peut-être corriger un peu certaines erreurs d'inclinaison, et re-situer les coordonnées spécifiques aux données (x, y) pour qu'elles correspondent à un système de coordonnées connu. Le résultat est censé produire des coordonnées projetées. Cela signifie simplement qu'il existe une procédure mathématique prenant (longitude, latitude) et les transformant en coordonnées connues: c'est ce qu'on appelle la "projection". "Unprojecting" fait l'inverse; Donc, si vous savez quelle projection est nécessaire, vous l'appliquez aux coordonnées affinées transformées (x, y) pour obtenir la latitude et la longitude.

Soit dit en passant, les valeurs des constantes a, b, c, d sont données par les entrées 1, 2, 4 et 5 du tableau de géotransformation.

Vous pouvez utiliser les éléments suivants:

import gdal, numpy as n
def coord(file):
    padfTransform = file.GetGeoTransform()
    indices = n.indices(file.ReadAsArray().shape)
    xp = padfTransform[0] + indices[1]*padfTransform[1] + indices[1]*padfTransform[2]   
    yp = padfTransform[3] + indices[0]*padfTransform[4] + indices[0]*padfTransform[5]  
    return xp,yp
file = gdal.Open('Yourfile.tif') # A GeoTiff file
x,y = coord(file)

coord renverra la longitude (x) et la latitude (y) de tous les pixels. Gardez à l'esprit que les coordonnées sont du coin gauche d'un pixel