Bagaimana cara mendapatkan koordinat sudut raster menggunakan binding Python GDAL?

30

Apakah ada cara untuk mendapatkan koordinat sudut (dalam derajat lat / panjang) dari file raster menggunakan binding Python gdal?

Beberapa pencarian online telah meyakinkan saya bahwa tidak ada, jadi saya telah mengembangkan pekerjaan dengan mem-parsing keluaran gdalinfo, ini agak mendasar tapi saya pikir mungkin menghemat waktu bagi orang yang mungkin kurang nyaman dengan python. Ini juga hanya berfungsi jika gdalinfo berisi koordinat geografis bersama dengan koordinat sudut, yang saya tidak percaya selalu demikian.

Ini solusi saya, apakah ada yang punya solusi yang lebih baik?

gdalinfo pada raster yang sesuai menghasilkan sesuatu seperti ini di tengah-tengah output:

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)

Kode ini akan berfungsi pada file yang terlihat seperti gdalinfo. Saya percaya kadang-kadang koordinat akan dalam derajat dan desimal, bukan derajat, menit dan detik; itu seharusnya sepele untuk menyesuaikan kode untuk situasi itu.

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

Dan itu memberi saya:

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

Jawaban:

29

Berikut cara lain untuk melakukannya tanpa memanggil program eksternal.

Apa yang dilakukan adalah mendapatkan koordinat keempat sudut dari geotransform dan memproyeksikannya kembali ke lon / lat menggunakan 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)

Beberapa kode dari metageta proyek, osr.CoordinateTransformation ide dari jawaban ini

pengguna2856
sumber
Keren terima kasih. Dan itu berfungsi terlepas dari apakah garis yang sesuai ada di output gdalinfo.
EddyTheB
Saya pikir ini akan lebih baik dengan tgt_srs = src_srs.CloneGeogCS (). Raster awal saya adalah gambar Mars, jadi menggunakan EPSG (4326) tidak ideal, tetapi CloneGeogCS () tampaknya hanya berubah dari proyeksi ke koordinat geografis.
EddyTheB
Tentunya. Saya telah memperbarui jawaban untuk menggunakan CloneGeogCS. Namun, saya hanya mencoba mendemonstrasikan penggunaan fungsi GetExtent dan ReprojectCoords. Anda dapat menggunakan apa pun yang Anda inginkan sebagai target srs.
user2856
Ya, terima kasih, saya belum pernah menemukan yang sebaliknya.
EddyTheB
Bagaimana jika Anda memiliki dataset yang tidak memiliki proyeksi dan menentukan RPC? Impor dari fungsi wkt gagal. Dimungkinkan untuk menghitung sejauh mana menggunakan transformator tapi saya bertanya-tanya apakah ada cara dengan metode di atas?
Thomas
41

Ini dapat dilakukan dalam baris kode yang jauh lebih sedikit

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, ulyadalah sudut kiri atas lrx,, lryadalah sudut kanan bawah

Pustaka osr (bagian dari gdal) dapat digunakan untuk mengubah titik ke sistem koordinat apa pun. Untuk satu titik:

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)

Untuk reproject citra raster seluruh akan menjadi masalah yang jauh lebih rumit, tapi GDAL> = 2.0 menawarkan solusi mudah untuk ini juga: gdal.Warp.

James
sumber
Ini adalah jawaban Pythonic untuk tingkat - solusi yang sama-Pythonic untuk proyeksi ulang itu akan luar biasa, Yang mengatakan - Saya menggunakan hasil di PostGIS, jadi saya hanya melewati batas yang tidak berubah dan ST_Transform(ST_SetSRID(ST_MakeBox2D([hasilnya] ),28355),4283)itu. (Satu berdalih - 'T' di src.GetGeoTransform()harus dikapitalisasi).
GT.
@ GT. Menambahkan sebuah contoh
James
1

Saya telah melakukan ini ... ini sedikit sulit dikodekan tetapi jika tidak ada perubahan dengan gdalinfo, ini akan berfungsi untuk gambar yang diproyeksikan 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
sumber
2
Ini cukup rapuh karena bergantung pada keduanya gdalinfotersedia di jalur pengguna (tidak selalu halnya, terutama pada windows) dan parsing output yang dicetak yang tidak memiliki antarmuka yang ketat - yaitu mengandalkan 'angka ajaib' untuk jarak yang benar. Tidak perlu ketika gdal menyediakan binding python komprehensif yang dapat melakukan ini dengan cara yang lebih eksplisit dan kuat
James
1

Jika raster Anda diputar, maka matematika menjadi sedikit lebih rumit, karena Anda perlu mempertimbangkan masing-masing dari enam koefisien transformasi affine. Pertimbangkan menggunakan affine untuk mengubah transformasi affine / geotransform yang diputar:

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|

Sekarang Anda dapat menghasilkan empat pasangan koordinat sudut:

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

Dan jika Anda juga membutuhkan batasan berbasis grid (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
sumber
0

Saya percaya pada versi yang lebih baru dari modul OSGEO / GDAL untuk python seseorang dapat langsung memanggil utilitas GDAL dari kode tanpa melibatkan panggilan sistem. misalnya alih-alih menggunakan subproses untuk menelepon:

gdalinfo orang dapat memanggil gdal.Info (the_name_of_the_file) untuk memiliki eksposur metadata file / anotasi

atau alih-alih menggunakan subproses untuk memanggil: gdalwarp seseorang dapat memanggil gdal. Lewati untuk melakukan warping pada raster.

Daftar utilitas GDAL saat ini tersedia sebagai fungsi internal: http://erouault.blogspot.com/2015/10/gdal-and-ogr-utilities-as-library.html

Sinooshka
sumber