对大疆御3E拍摄照片赋予坐标系并旋转

问题描述

进行植被覆盖度验证时，需采集验证点的植被覆盖情况，但无人机拍摄的照片缺少坐标系，无法进行对比验证。

解决方案

赋予照片坐标系

在gdal中对影像赋予坐标系主要参数为仿射六参数：左上角坐标，xy分辨率和旋转信息，故只需解决该6个参数信息即可将照片初步赋予坐标系，本次以赋予WGS-84坐标系为例。

左上角坐标计算

查看照片属性中的详细信息，其中记录了照片拍摄时的GPS信息（中心坐标）或对照片使用记事本打开，其中 drone-dji:GpsLatitude="“和drone-dji:GpsLongitude=”"两个参数也记录了经纬度信息。根据该GPS信息进行左上角坐标计算，方法为将该坐标视为图像中心坐标，无人机拍摄的照片分辨率已知（3cm），根据中心坐标的行列号与左上角坐标的行列号、xy分辨率计算出左上角坐标计算。

import re
from osgeo import gdal
import os

def png_deal(tif_src, png_file, out_folder, out_name):
    '''
    :param tif_src:正射影像的dataset，用于获取xy分辨率（如果这里没有正射影像，可自行设定分辨率）
    :param png_file:待校正照片的绝对路径
    :param out_folder:输出文件夹
    :param out_name:输出名字
    :return:
    '''
    tif_geo = tif_src.GetGeoTransform()
    tif_x_res = tif_geo[1]
    tif_y_res = tif_geo[5]
    png_src = gdal.Open(png_file)
    # get png metadata
    png_meta = png_src.GetMetadata()
    # get latitude longitude x y
    png_lon = png_meta['EXIF_GPSLongitude']
    png_lat = png_meta['EXIF_GPSLatitude']
    png_x = png_src.RasterXSize
    png_y = png_src.RasterYSize


    # change mile to degree
    png_x_res = tif_x_res / (2 * math.pi * 6371004) * 360;
    png_y_res = tif_y_res / (2 * math.pi * 6371004) * 360;

    lon_str = re.findall(r'\d+', png_lon)
    lon_decimal = float(lon_str[0]) + float(lon_str[1])/60 + (float(lon_str[2]) )/3600
    lat_str = re.findall(r'\d+', png_lat)
    lat_decimal = float(lat_str[0]) + float(lat_str[1])/60 + (float(lat_str[2]))/3600

    # caculate the left-top coordinate
    png_left_coord = [lon_decimal - png_x_res * png_x/2, lat_decimal - png_y_res * png_y/2]
    out_png_geom = [png_left_coord[0], png_x_res, 0.0, png_left_coord[1], 0.0, png_y_res]

    srs = osr.SpatialReference()
    srs.ImportFromEPSG(4326)

    # export the png tif
    driver = gdal.GetDriverByName('GTiff')
    # out_tif = driver.Create(name=out_name, ysize=tar_y, xsize=tar_x, bands=tar_bandnum, eType=tar_datatype)
    out_tif = driver.Create(out_name, png_x, png_y, png_src.RasterCount, eType=png_src.GetRasterBand(1).DataType)
    for i in range(png_src.RasterCount):
        data = png_src.GetRasterBand(i+1).ReadAsArray()
        band = out_tif.GetRasterBand(i+1).WriteArray(data)
        del data, band

    out_tif.SetProjection(srs.ExportToWkt())
    out_tif.SetGeoTransform(out_png_geom)
    out_tif.FlushCache()

    return out_tif

对照片进行旋转

对照片使用记事本打开，其中drone-dji:FlightYawDegree=""记录了照片的旋转信息，这里借助Affrine库实现对影像的旋转

# 该函数主要获取照片的偏航角
def Get_Image_Yaw_angle(file_path):
    """
    :param file_path: 输入图片路径
    :return: 图片的偏航角
    """
    # 获取图片偏航角
    b = b"\x3c\x2f\x72\x64\x66\x3a\x44\x65\x73\x63\x72\x69\x70\x74\x69\x6f\x6e\x3e"
    a = b"\x3c\x72\x64\x66\x3a\x44\x65\x73\x63\x72\x69\x70\x74\x69\x6f\x6e\x20"
    img = open(file_path, 'rb')
    data = bytearray()
    dj_data_dict = {}
    flag = False
    for line in img.readlines():
        if a in line:
            flag = True
        if flag:
            data += line
        if b in line:
            break
    if len(data) > 0:
        data = str(data.decode('ascii'))
        lines = list(filter(lambda x: 'drone-dji:' in x, data.split("\n")))
        for d in lines:
            d = d.strip()[10:]
            key, value = d.split("=")
            dj_data_dict[key] = value
    # print("Image_yaw",dj_data_dict["FlightYawDegree"][1:-1])
    return float(dj_data_dict["FlightYawDegree"][1:-1])

# 获取中心像元行列号
def raster_center(ds):
    """This function return the pixel coordinates of the raster center
    """

    # We get the size (in pixels) of the raster
    # using gdal
    width, height = ds.RasterXSize, ds.RasterYSize

    # We calculate the middle of raster
    xmed = width / 2
    ymed = height / 2
    tar_geom = ds.GetGeoTransform()
    tar_Xp = tar_geom [0] + xmed * tar_geom[1] + ymed * tar_geom[2]
    tar_Yp = tar_geom [3] + xmed * tar_geom[4] + ymed * tar_geom[5]

    return [tar_Xp, tar_Yp]
# 对照片进行旋转
def rotate_gt(affine_matrix, angle, pivot=None):
    """This function generate a rotated affine matrix
    """

    affine_src = Affine.from_gdal(*affine_matrix)
    # We made the rotation. For this we calculate a rotation matrix,
    # with the rotation method and we combine it with the original affine matrix
    # Be carful, the star operator (*) is surcharged by Affine package. He make
    # a matrix multiplication, not a basic multiplication
    affine_dst = affine_src * affine_src.rotation(angle, pivot)
    # We retrun the rotated matrix in gdal format
    return affine_dst.to_gdal()

yaw_angle = Get_Image_Yaw_angle(jpg_file)
after_coordinate_file = r'输入赋予坐标系后的照片路径'
dataset_src = gdal.Open(after_coordinate_file)
out_fin_folder = r'输出文件夹'
out_file = r'输出绝对路径'
# 创建输出路径
driver = gdal.GetDriverByName('GTiff')
datase_dst = driver.CreateCopy(out_file, dataset_src, strict=0)
gt_affine = dataset_src.GetGeoTransform()
center = raster_center(dataset_src)
# 进行旋转输出
datase_dst.SetGeoTransform(rotate_gt(gt_affine, yaw_angle, center))

结果如下
原始照片：

校正后结果：

参考：
https://blog.csdn.net/m0_56729804/article/details/131695618