Raster georeferensi menggunakan GDAL dan Python?

10

Saya ingin melakukan georeferensi raster menggunakan pythondan GDAL. Pendekatan saya saat ini adalah menelepon gdal_translatedan gdalwarpmenggunakan os.systemdan daftar titik kontrol tanah yang jelek. Saya benar-benar ingin cara untuk melakukan ini secara asli di dalam python.

Ini adalah proses saat ini yang saya gunakan:

import os

os.system('gdal_translate -of GTiff -gcp 1251.92 414.538 -7.9164e+06 5.21094e+06 -gcp 865.827 107.699 -7.91651e+06 5.21104e+06  "inraster.tif" "outraster1.tif"')
os.system('gdalwarp -r bilinear -tps -co COMPRESS=NONE "outraster2.tif" "outraster3.tif"')

Ada pertanyaan dan jawaban sebelumnya dari 2012 yang menyatakan bahwa gdal_translatedapat diakses setelah mengimpor gdal. Saya tidak yakin apakah sudah usang, atau apakah itu salah tetapi ketika saya menjalankan from osgeo import gdalsaya tidak melihat gdal.gdal_translatesebagai opsi.

Saya tidak tahu apakah itu ada tetapi saya akan senang jika saya bisa menerjemahkan dan memproyeksikan ulang raster dengan cara pythonic. Sebagai contoh:

# translate
gcp_points = [(1251.92, 414.538), (-7.9164e+06, 5.21094e+06)]
gdal.gdal_translate(in_raster, gcp_points, out_raster1)

# warp
gdal.gdalwarp(out_raster1, out_raster2, 'bilinear', args*)

Apakah pendekatan seperti itu mungkin?

djq
sumber
1
gdal_translate dan gdalwarp dan utilitas gdal lainnya bukanlah paket / modul python, mereka adalah executable mandiri. Anda dapat menggunakan sistem os.s atau (lebih disukai) proses. Buka untuk memanggil mereka.
user2856
@ Lukas jadi apakah tidak ada cara untuk berinteraksi dengan utilitas gdal ini? Saya akan senang menjelajahi cara menulis pembungkus python untuk ini jika ada yang tidak ada.
djq
Anda dapat menggunakan API python gdal untuk melakukan hampir semua yang dapat dilakukan utilitas gdal_translate dan gdalwarp, itu hanya sedikit lebih rumit. Lihat gdal. AutoCreateWarpedVRT dan gdal Driver.CreateCopy.
user2856

Jawaban:

17

Berikut adalah contoh yang kira-kira sesuai dengan apa yang Anda minta. Parameter utama adalah geotransformlarik yang gdal gunakan untuk menggambarkan lokasi raster (posisi, skala piksel, dan kemiringan) dan epsgkode proyeksi. Dengan itu, kode berikut harus melakukan georeferensi raster dengan benar dan menentukan proyeksi.

Saya tidak menguji sebanyak ini, tetapi tampaknya berhasil bagi saya. Anda harus memastikan koordinat yang saya masukkan benar dan mungkin mengubah proyeksi dan skala pixel / ukuran. Semoga ini bisa membantu.

from osgeo import gdal, osr

src_filename ='/path/to/source.tif'
dst_filename = '/path/to/destination.tif'

# Opens source dataset
src_ds = gdal.Open(src_filename)
format = "GTiff"
driver = gdal.GetDriverByName(format)

# Open destination dataset
dst_ds = driver.CreateCopy(dst_filename, src_ds, 0)

# Specify raster location through geotransform array
# (uperleftx, scalex, skewx, uperlefty, skewy, scaley)
# Scale = size of one pixel in units of raster projection
# this example below assumes 100x100
gt = [-7916400, 100, 0, 5210940, 0, -100]

# Set location
dst_ds.SetGeoTransform(gt)

# Get raster projection
epsg = 3857
srs = osr.SpatialReference()
srs.ImportFromEPSG(epsg)
dest_wkt = srs.ExportToWkt()

# Set projection
dst_ds.SetProjection(dest_wkt)

# Close files
dst_ds = None
src_ds = None
yellowcap
sumber
2
Itu benar tetapi jika Anda ingin lebih akurat, dalam hal rotasi, Anda harus menggunakan Affine Transformation dan menghitung parameter yang tepat dengan paket Affine ( unduh dari sini ). Penggunaan: 'Affine.scale (PARAMETER) * Affine.rotation (ANGLE)'. PARAMETER berasal dari rasio piksel dari dua gambar, Anda dapat menggunakan 'gdalinfo' atau PIL untuk mengetahui ukuran piksel, ANGLE adalah sudutnya.
CaMa