Как установить географическую привязку аэрофотоснимка без привязки, используя наземные контрольные точки в Python - PullRequest
Новая
       25

Как установить географическую привязку аэрофотоснимка без привязки, используя наземные контрольные точки в Python

2 голосов
/ 15 апреля 2019

У меня есть серия аэрофотоснимков без ссылок, которые я бы хотел использовать для привязки с помощью python.Изображения пространственно идентичны (на самом деле они представляют собой кадры, извлеченные из видео), и я получил для них контрольные точки с помощью географической привязки одного кадра в ArcMap вручную.Я хотел бы применить полученные наземные контрольные точки ко всем последующим изображениям и в результате получить гео-формат или файл JPEG с соответствующим файлом мира (.jgw) для каждого обработанного изображения.Я знаю, что это возможно сделать с помощью arcpy, но у меня нет доступа к arcpy, и я действительно хотел бы использовать бесплатный модуль с открытым исходным кодом, если это возможно.

Моя система координат - NZGD2000 (epsg 2193), ивот таблица контрольных точек, которые я хочу применить к своим изображениям:

176.412984, -310.977264, 1681255.524654, 6120217.357425

160.386905, -141.487145, 1681158.424227, 6120406.821249 * 100 *, -310.547238, 1681556.948690, 6120335.658359

Вот пример изображения: https://imgur.com/a/9ThHtOz

Я прочитал много информации о GDAL и растерио, но у меня нет опыта работы сих, и я не могу адаптировать биты кода, которые я нашел в моей конкретной ситуации.

Попытка Растерио: 

import cv2
from rasterio.warp import reproject
from rasterio.control import GroundControlPoint
from fiona.crs import from_epsg

img = cv2.imread("Example_image.jpg")

# Creating ground control points (not sure if I got the order of variables right):
points = [(GroundControlPoint(176.412984, -310.977264, 1681255.524654, 6120217.357425)),
          (GroundControlPoint(160.386905, -141.487145, 1681158.424227, 6120406.821253)),
          (GroundControlPoint(433.204947, -310.547238, 1681556.948690, 6120335.658359))]

# The function requires a parameter "destination", but I'm not sure what to put there.
#   I'm guessing this may not be the right function to use
reproject(img, destination, src_transform=None, gcps=points, src_crs=from_epsg(2193),
                        src_nodata=None, dst_transform=None, dst_crs=from_epsg(2193), dst_nodata=None,
                        src_alpha=0, dst_alpha=0, init_dest_nodata=True, warp_mem_limit=0)

Попытка GDAL: 

from osgeo import gdal 
import osr

inputImage = "Example_image.jpg"
outputImage = "image_gdal.jpg"

dataset = gdal.Open(inputImage) 
I = dataset.ReadAsArray(0,0,dataset.RasterXSize,dataset.RasterYSize)

outdataset = gdal.GetDriverByName('GTiff') 
output_SRS = osr.SpatialReference() 
output_SRS.ImportFromEPSG(2193) 
outdataset = outdataset.Create(outputImage,dataset.RasterXSize,dataset.RasterYSize,I.shape[0]) 
for nb_band in range(I.shape[0]):
    outdataset.GetRasterBand(nb_band+1).WriteArray(I[nb_band,:,:])

# Creating ground control points (not sure if I got the order of variables right):
gcp_list = [] 
gcp_list.append(gdal.GCP(176.412984, -310.977264, 1681255.524654, 6120217.357425))
gcp_list.append(gdal.GCP(160.386905, -141.487145, 1681158.424227, 6120406.821253))
gcp_list.append(gdal.GCP(433.204947, -310.547238, 1681556.948690, 6120335.658359))

outdataset.SetProjection(srs.ExportToWkt()) 
wkt = outdataset.GetProjection() 
outdataset.SetGCPs(gcp_list,wkt)

outdataset = None

Не знаюЯ не знаю, как заставить работать приведенный выше код, и я был бы очень признателен за любую помощь в этом.

1 Ответ

1 голос
/ 17 апреля 2019

Для вашего метода gdal должно работать только использование gdal.Warp с набором выходных данных, например, 

outdataset.SetProjection(srs.ExportToWkt()) 
wkt = outdataset.GetProjection() 
outdataset.SetGCPs(gcp_list,wkt)
gdal.Warp("output_name.tif", outdataset, dstSRS='EPSG:2193', format='gtiff')

Это создаст новый файл output_name.tif.

...