Comment obtenir des coordonnées de coin raster à l'aide de liaisons Python GDAL?

30

Existe-t-il un moyen d'obtenir les coordonnées des coins (en degrés lat / long) à partir d'un fichier raster en utilisant les liaisons Python de gdal?

Quelques recherches en ligne m'ont convaincu qu'il n'y en avait pas, j'ai donc développé un moyen de contourner la sortie de gdalinfo, c'est un peu basique mais je pensais que cela pourrait faire gagner du temps aux personnes qui seraient moins à l'aise avec python. Cela ne fonctionne également que si gdalinfo contient les coordonnées géographiques ainsi que les coordonnées des coins, ce qui, je pense, n'est pas toujours le cas.

Voici ma solution, quelqu'un a-t-il de meilleures solutions?

gdalinfo sur un raster approprié se traduit par quelque chose comme ceci à mi-chemin de la sortie:

Corner Coordinates:
Upper Left  (  -18449.521, -256913.934) (137d 7'21.93"E,  4d20'3.46"S)
Lower Left  (  -18449.521, -345509.683) (137d 7'19.32"E,  5d49'44.25"S)
Upper Right (   18407.241, -256913.934) (137d44'46.82"E,  4d20'3.46"S)
Lower Right (   18407.241, -345509.683) (137d44'49.42"E,  5d49'44.25"S)
Center      (     -21.140, -301211.809) (137d26'4.37"E,  5d 4'53.85"S)

Ce code fonctionnera sur les fichiers qui ressemblent à gdalinfo. Je crois que parfois les coordonnées seront en degrés et décimales, plutôt qu'en degrés, minutes et secondes; il devrait être trivial d'ajuster le code pour cette situation.

import numpy as np
import subprocess

def GetCornerCoordinates(FileName):
    GdalInfo = subprocess.check_output('gdalinfo {}'.format(FileName), shell=True)
    GdalInfo = GdalInfo.split('/n') # Creates a line by line list.
    CornerLats, CornerLons = np.zeros(5), np.zeros(5)
    GotUL, GotUR, GotLL, GotLR, GotC = False, False, False, False, False
    for line in GdalInfo:
        if line[:10] == 'Upper Left':
            CornerLats[0], CornerLons[0] = GetLatLon(line)
            GotUL = True
        if line[:10] == 'Lower Left':
            CornerLats[1], CornerLons[1] = GetLatLon(line)
            GotLL = True
        if line[:11] == 'Upper Right':
            CornerLats[2], CornerLons[2] = GetLatLon(line)
            GotUR = True
        if line[:11] == 'Lower Right':
            CornerLats[3], CornerLons[3] = GetLatLon(line)
            GotLR = True
        if line[:6] == 'Center':
            CornerLats[4], CornerLons[4] = GetLatLon(line)
            GotC = True
        if GotUL and GotUR and GotLL and GotLR and GotC:
            break
    return CornerLats, CornerLons 

def GetLatLon(line):
    coords = line.split(') (')[1]
    coords = coords[:-1]
    LonStr, LatStr = coords.split(',')
    # Longitude
    LonStr = LonStr.split('d')    # Get the degrees, and the rest
    LonD = int(LonStr[0])
    LonStr = LonStr[1].split('\'')# Get the arc-m, and the rest
    LonM = int(LonStr[0])
    LonStr = LonStr[1].split('"') # Get the arc-s, and the rest
    LonS = float(LonStr[0])
    Lon = LonD + LonM/60. + LonS/3600.
    if LonStr[1] in ['W', 'w']:
        Lon = -1*Lon
    # Same for Latitude
    LatStr = LatStr.split('d')
    LatD = int(LatStr[0])
    LatStr = LatStr[1].split('\'')
    LatM = int(LatStr[0])
    LatStr = LatStr[1].split('"')
    LatS = float(LatStr[0])
    Lat = LatD + LatM/60. + LatS/3600.
    if LatStr[1] in ['S', 's']:
        Lat = -1*Lat
    return Lat, Lon

FileName = Image.cub
# Mine's an ISIS3 cube file.
CornerLats, CornerLons = GetCornerCoordinates(FileName)
# UpperLeft, LowerLeft, UpperRight, LowerRight, Centre
print CornerLats
print CornerLons

Et cela me donne:

[-4.33429444 -5.82895833 -4.33429444 -5.82895833 -5.081625  ] 
[ 137.12275833  137.12203333  137.74633889  137.74706111  137.43454722]
EddyTheB
la source

Réponses:

29

Voici une autre façon de le faire sans appeler un programme externe.

Cela permet d'obtenir les coordonnées des quatre coins de la géotransformation et de les reprojeter en lon / lat en utilisant osr.CoordinateTransformation.

from osgeo import gdal,ogr,osr

def GetExtent(gt,cols,rows):
    ''' Return list of corner coordinates from a geotransform

        @type gt:   C{tuple/list}
        @param gt: geotransform
        @type cols:   C{int}
        @param cols: number of columns in the dataset
        @type rows:   C{int}
        @param rows: number of rows in the dataset
        @rtype:    C{[float,...,float]}
        @return:   coordinates of each corner
    '''
    ext=[]
    xarr=[0,cols]
    yarr=[0,rows]

    for px in xarr:
        for py in yarr:
            x=gt[0]+(px*gt[1])+(py*gt[2])
            y=gt[3]+(px*gt[4])+(py*gt[5])
            ext.append([x,y])
            print x,y
        yarr.reverse()
    return ext

def ReprojectCoords(coords,src_srs,tgt_srs):
    ''' Reproject a list of x,y coordinates.

        @type geom:     C{tuple/list}
        @param geom:    List of [[x,y],...[x,y]] coordinates
        @type src_srs:  C{osr.SpatialReference}
        @param src_srs: OSR SpatialReference object
        @type tgt_srs:  C{osr.SpatialReference}
        @param tgt_srs: OSR SpatialReference object
        @rtype:         C{tuple/list}
        @return:        List of transformed [[x,y],...[x,y]] coordinates
    '''
    trans_coords=[]
    transform = osr.CoordinateTransformation( src_srs, tgt_srs)
    for x,y in coords:
        x,y,z = transform.TransformPoint(x,y)
        trans_coords.append([x,y])
    return trans_coords

raster=r'somerasterfile.tif'
ds=gdal.Open(raster)

gt=ds.GetGeoTransform()
cols = ds.RasterXSize
rows = ds.RasterYSize
ext=GetExtent(gt,cols,rows)

src_srs=osr.SpatialReference()
src_srs.ImportFromWkt(ds.GetProjection())
#tgt_srs=osr.SpatialReference()
#tgt_srs.ImportFromEPSG(4326)
tgt_srs = src_srs.CloneGeogCS()

geo_ext=ReprojectCoords(ext,src_srs,tgt_srs)

Du code du projet metageta, idée de transformation osr.CoordinateTransformation de cette réponse

user2856
la source
Super merci. Et cela fonctionne indépendamment du fait que les lignes appropriées existent dans la sortie gdalinfo.
EddyTheB
Je pense que ce sera mieux avec tgt_srs = src_srs.CloneGeogCS (). Mes rasters initiaux sont des images de Mars, donc l'utilisation d'EPSG (4326) n'est pas idéale, mais CloneGeogCS () semble simplement passer des coordonnées projetées aux coordonnées géographiques.
EddyTheB
Pour sûr. J'ai mis à jour la réponse pour utiliser CloneGeogCS. Cependant, j'essayais simplement de démontrer l'utilisation des fonctions GetExtent et ReprojectCoords. Vous pouvez utiliser tout ce que vous voulez comme cible srs.
user2856
Oui, merci, je n'en aurais jamais trouvé autrement.
EddyTheB
Que se passe-t-il si vous avez un ensemble de données qui n'a pas de projection et spécifie les RPC? L'importation à partir de la fonction wkt échoue. Il est possible de calculer l'étendue à l'aide d'un transformateur mais je me demandais s'il y avait un moyen avec la méthode ci-dessus?
Thomas
41

Cela peut être fait en beaucoup moins de lignes de code

src = gdal.Open(path goes here)
ulx, xres, xskew, uly, yskew, yres  = src.GetGeoTransform()
lrx = ulx + (src.RasterXSize * xres)
lry = uly + (src.RasterYSize * yres)

ulx, ulyEst le coin supérieur gauche, lrx, lryest le coin inférieur droit

La bibliothèque osr (partie de gdal) peut être utilisée pour transformer les points en n'importe quel système de coordonnées. Pour un seul point:

from osgeo import ogr
from osgeo import osr

# Setup the source projection - you can also import from epsg, proj4...
source = osr.SpatialReference()
source.ImportFromWkt(src.GetProjection())

# The target projection
target = osr.SpatialReference()
target.ImportFromEPSG(4326)

# Create the transform - this can be used repeatedly
transform = osr.CoordinateTransformation(source, target)

# Transform the point. You can also create an ogr geometry and use the more generic `point.Transform()`
transform.TransformPoint(ulx, uly)

Pour reprojeter une image entière de trame serait une question beaucoup plus compliquée, mais GDAL> = 2.0 offre une solution simple pour cela aussi: gdal.Warp.

James
la source
C'est la réponse Pythonic pour l'étendue - une solution tout aussi Pythonic pour la reprojection aurait été géniale, Cela dit - j'utilise les résultats dans PostGIS, donc je passe juste l'étendue non transformée et ST_Transform(ST_SetSRID(ST_MakeBox2D([les résultats] ),28355),4283). (Un petit problème - le «T» src.GetGeoTransform()doit être en majuscule).
GT.
@GT. Ajout d'un exemple
James
1

Je l'ai fait de cette façon ... c'est un peu codé en dur mais si rien ne change avec gdalinfo, cela fonctionnera pour les images projetées UTM!

imagefile= /pathto/image
p= subprocess.Popen(["gdalinfo", "%s"%imagefile], stdout=subprocess.PIPE)
out,err= p.communicate()
ul= out[out.find("Upper Left")+15:out.find("Upper Left")+38]
lr= out[out.find("Lower Right")+15:out.find("Lower Right")+38]
Emiliano
la source
2
Ceci est assez fragile car il dépend à la fois d' gdalinfoêtre disponible sur le chemin d'un utilisateur (pas toujours le cas, en particulier sur Windows) et d'analyser une sortie imprimée qui n'a pas d'interface stricte - c'est-à-dire de s'appuyer sur ces `` nombres magiques '' pour un espacement correct. Il n'est pas nécessaire lorsque gdal fournit des liaisons python complètes qui peuvent le faire de manière plus explicite et robuste
James
1

Si votre raster est pivoté, les calculs deviennent un peu plus compliqués, car vous devez considérer chacun des six coefficients de transformation affine. Envisagez d'utiliser l' affine pour transformer une transformation / géotransformation affine tournée:

from affine import Affine

# E.g., from a GDAL DataSet object:
# gt = ds.GetGeoTransform()
# ncol = ds.RasterXSize
# nrow = ds.RasterYSize

# or to work with a minimal example
gt = (100.0, 17.320508075688775, 5.0, 200.0, 10.0, -8.660254037844387)
ncol = 10
nrow = 15

transform = Affine.from_gdal(*gt)
print(transform)
# | 17.32, 5.00, 100.00|
# | 10.00,-8.66, 200.00|
# | 0.00, 0.00, 1.00|

Vous pouvez maintenant générer les quatre paires de coordonnées de coin:

c0x, c0y = transform.c, transform.f  # upper left
c1x, c1y = transform * (0, nrow)     # lower left
c2x, c2y = transform * (ncol, nrow)  # lower right
c3x, c3y = transform * (ncol, 0)     # upper right

Et si vous avez également besoin des limites basées sur la grille (xmin, ymin, xmax, ymax):

xs = (c0x, c1x, c2x, c3x)
ys = (c0y, c1y, c2y, c3y)
bounds = min(xs), min(ys), max(xs), max(ys)
Mike T
la source
0

Je crois que dans les versions plus récentes du module OSGEO / GDAL pour python, on peut appeler directement les utilitaires GDAL à partir du code sans impliquer d'appels système. par exemple au lieu d'utiliser un sous-processus pour appeler:

gdalinfo on peut appeler gdal.Info (the_name_of_the_file) pour avoir une exposition des métadonnées / annotations du fichier

ou au lieu d'utiliser un sous-processus pour appeler: gdalwarp on peut appeler gdal.Warp pour effectuer la déformation sur un raster.

La liste des utilitaires GDAL actuellement disponibles en tant que fonction interne: http://erouault.blogspot.com/2015/10/gdal-and-ogr-utilities-as-library.html

Sinooshka
la source