如何使用python中的地面控制点对未参考的航拍图像进行地理配准答案

【问题标题】：How to georeference an unreferenced aerial image using ground control points in python如何使用python中的地面控制点对未参考的航拍图像进行地理配准
【发布时间】：2019-04-15 02:47:00
【问题描述】：

我有一系列未参考的航拍图像，我想使用 python 进行地理配准。这些图像在空间上是相同的（它们实际上是从视频中提取的帧），我通过在 ArcMap 中手动对一帧进行地理配准获得了它们的地面控制点。我想将获得的地面控制点应用于所有后续图像，并因此为每个处理后的图像获得一个带有相应世界文件 (.jgw) 的 geo-tiff 或 jpeg 文件。我知道使用 arcpy 可以做到这一点，但我无法访问 arcpy，如果可能的话，我真的很想使用免费的开源模块。

我的坐标系是 NZGD2000 (epsg 2193)，这里是我希望应用于我的图像的控制点表：

176.412984、-310.977264、1681255.524654、6120217.357425

160.386905、-141.487145、1681158.424227、6120406.821253

433.204947、-310.547238、1681556.948690、6120335.658359

这是一个示例图片：https://imgur.com/a/9ThHtOz

我已经阅读了很多关于 GDAL 和 rasterio 的信息，但我对它们没有任何经验，并且无法根据我的特定情况调整我发现的代码。

光栅尝试：

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

我不太清楚如何使上述代码工作，如果能提供任何帮助，我将不胜感激。

【问题讨论】：

标签： python image gis gdal rasterio

【解决方案1】：

我最终阅读了一本书“使用 Python 进行地理处理”，并最终找到了适合我的解决方案。这是我适应我的问题的代码：

import shutil
from osgeo import gdal, osr

orig_fn = 'image.tif'
output_fn = 'output.tif'

# Create a copy of the original file and save it as the output filename:
shutil.copy(orig_fn, output_fn)
# Open the output file for writing for writing:
ds = gdal.Open(output_fn, gdal.GA_Update)
# Set spatial reference:
sr = osr.SpatialReference()
sr.ImportFromEPSG(2193) #2193 refers to the NZTM2000, but can use any desired projection

# Enter the GCPs
#   Format: [map x-coordinate(longitude)], [map y-coordinate (latitude)], [elevation],
#   [image column index(x)], [image row index (y)]
gcps = [gdal.GCP(1681255.524654, 6120217.357425, 0, 176.412984, 310.977264),
gdal.GCP(1681158.424227, 6120406.821253, 0, 160.386905, 141.487145),
gdal.GCP(1681556.948690, 6120335.658359, 0, 433.204947, 310.547238)]

# Apply the GCPs to the open output file:
ds.SetGCPs(gcps, sr.ExportToWkt())

# Close the output file in order to be able to work with it in other programs:
ds = None

【讨论】：

感谢您的回答！我有一个类似的问题，这有助于我了解如何解决它。我知道您可能会对此做出回应，但您知道是否可以仅使用一个 GCP 来执行此操作。我也知道像素的分辨率（以米为单位）。如果您有兴趣回答，这是我的问题：gis.stackexchange.com/questions/386863/….
@salRad 看起来你已经想通了，但如果你仍然想要我的意见，请告诉我:)
@Kat 感谢您的回答，您让我的生活更轻松。对于可能想要使用 jpg 而不是 tiff 的任何人，整个过程都可以完美运行，但是您需要在 gdal.Open 中取出 gdal.GA_Update 选项才能读取 jpg 图像！干杯

【解决方案2】：

对于您的 gdal 方法，只需将 gdal.Warp 与 outdataset 一起使用即可，例如

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。

【讨论】：

我最初挂断srs.ExportToWkt()，收到以下错误：“NameError: name 'srs' is not defined”。我刚刚将这段代码更改为output_SRS.ExportToWkt()，现在该部分似乎运行了。现在整个 GDAL 脚本运行并生成一个“image_gdal.jpg”文件，该文件没有地理参考，也没有任何与之关联的世界文件。我按照您的建议添加了gdal.Warp("output_name.tif", outdataset, dstSRS='EPSG:2193', format='gtiff')，但找不到“output_name.tif” - 它似乎没有生成。

【解决方案3】：

作为@Kat回答的补充，为了避免原始图像文件的质量损失并将nodata-value设置为0，可以使用以下内容。

#Load the original file
src_ds = gdal.Open(orig_fn)

#Create tmp dataset saved in memory
driver = gdal.GetDriverByName('MEM')
tmp_ds = driver.CreateCopy('', src_ds, strict=0)

#
# ... setting GCP....
#

# Setting no data for all bands
for i in range(1, tmp_ds.RasterCount + 1):
    f = tmp_ds.GetRasterBand(i).SetNoDataValue(0)

# Saving as file
driver = gdal.GetDriverByName('GTiff')
ds = driver.CreateCopy(output_fn, tmp_ds, strict=0)

【讨论】：