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

本文主要是介绍python学习——对大疆御3E拍摄照片赋予坐标系并旋转，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！

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

问题描述

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

解决方案

赋予照片坐标系

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

左上角坐标计算

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

import re
from osgeo import gdal
import osdef 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 metadatapng_meta = png_src.GetMetadata()# get latitude longitude x ypng_lon = png_meta['EXIF_GPSLongitude']png_lat = png_meta['EXIF_GPSLatitude']png_x = png_src.RasterXSizepng_y = png_src.RasterYSize# change mile to degreepng_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]) )/3600lat_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 coordinatepng_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 tifdriver = 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, bandout_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 = Falsefor line in img.readlines():if a in line:flag = Trueif flag:data += lineif b in line:breakif 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 gdalwidth, height = ds.RasterXSize, ds.RasterYSize# We calculate the middle of rasterxmed = width / 2ymed = height / 2tar_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 multiplicationaffine_dst = affine_src * affine_src.rotation(angle, pivot)# We retrun the rotated matrix in gdal formatreturn 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))