我有一系列未参考的航空影像,希望使用Python进行地理参考。这些图像在空间上是相同的(实际上它们是从视频中提取的帧),我通过手动在ArcMap中对其中一个帧进行地理参考来获取了它们的地面控制点。我想将我获得的地面控制点应用于所有后续图像,并因此获得一个geo-tiff或每个处理图像的相应世界文件(.jgw)的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的信息,但我没有任何使用它们的经验,并且无法将我找到的代码片段适应于我的特定情况。
尝试使用Rasterio:
我的坐标系统是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的信息,但我没有任何使用它们的经验,并且无法将我找到的代码片段适应于我的特定情况。
尝试使用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
我不是很清楚如何使上面的代码工作,非常感谢您对此提供任何帮助。