Python实战：用GDAL给大疆御3E照片添加WGS-84坐标系（附完整代码）

Python实战：用GDAL给大疆御3E照片添加WGS-84坐标系（附完整代码）

无人机航拍影像在地理信息系统（GIS）分析中扮演着越来越重要的角色。大疆御3E作为行业级无人机，其拍摄的高分辨率照片常被用于精准农业、环境监测等领域。然而，这些照片往往缺乏地理坐标系信息，导致无法直接用于空间分析。本文将手把手教你如何用Python和GDAL库解决这个问题。

1. 准备工作与环境配置

处理无人机影像需要特定的工具链和库支持。首先确保你的Python环境已安装以下关键组件：

• GDAL：地理数据抽象库，支持多种栅格和矢量格式
• Pillow：图像处理基础库
• NumPy：科学计算基础包

推荐使用conda安装GDAL，可以避免复杂的依赖问题：

conda install -c conda-forge gdal

对于Windows用户，如果遇到安装问题，可以考虑从GIS Internals下载预编译的GDAL二进制包。安装完成后，验证GDAL版本：

from osgeo import gdal print(gdal.__version__)

提示：大疆御3E拍摄的照片通常包含EXIF元数据，其中就包括关键的GPS位置信息。这些信息将以特定格式嵌入在JPEG文件中。

2. 解析无人机照片元数据

大疆无人机拍摄的照片会在EXIF中存储丰富的元数据，我们需要从中提取出关键的定位信息。以下是典型的元数据结构：

使用Python提取这些信息的核心代码如下：

from osgeo import gdal import re def extract_metadata(image_path): dataset = gdal.Open(image_path) if dataset is None: raise ValueError("无法打开图像文件") metadata = dataset.GetMetadata() lat = metadata.get('EXIF_GPSLatitude', '') lon = metadata.get('EXIF_GPSLongitude', '') # 转换度分秒格式为十进制 def dms_to_decimal(dms): parts = re.findall(r'\d+', dms) if len(parts) >= 3: return float(parts[0]) + float(parts[1])/60 + float(parts[2])/3600 return 0.0 return { 'latitude': dms_to_decimal(lat), 'longitude': dms_to_decimal(lon), 'width': dataset.RasterXSize, 'height': dataset.RasterYSize }

3. 构建地理坐标系

WGS-84（World Geodetic System 1984）是全球通用的地理坐标系。要为无人机照片添加这种坐标系，需要计算六个关键参数：

1. 左上角X坐标（经度）
2. X方向分辨率
3. 旋转参数（通常为0）
4. 左上角Y坐标（纬度）
5. 旋转参数（通常为0）
6. Y方向分辨率

计算这些参数的数学原理如下：

左上角经度 = 中心经度 - (图像宽度/2) * 经度分辨率 左上角纬度 = 中心纬度 - (图像高度/2) * 纬度分辨率

以下是完整的坐标转换函数：

from osgeo import osr def create_geotransform(metadata, resolution=0.00001): """ 根据元数据创建仿射变换参数 :param metadata: 包含经纬度和图像尺寸的字典 :param resolution: 每个像素代表的经纬度度数 :return: 六参数仿射变换数组 """ center_lon = metadata['longitude'] center_lat = metadata['latitude'] # 计算左上角坐标 top_left_lon = center_lon - (metadata['width']/2) * resolution top_left_lat = center_lat + (metadata['height']/2) * resolution # 注意纬度分辨率为负 return [top_left_lon, resolution, 0, top_left_lat, 0, -resolution] def assign_coordinate_system(input_path, output_path): metadata = extract_metadata(input_path) geotransform = create_geotransform(metadata) src_ds = gdal.Open(input_path) driver = gdal.GetDriverByName('GTiff') # 创建输出文件 dst_ds = driver.CreateCopy(output_path, src_ds) dst_ds.SetGeoTransform(geotransform) # 设置WGS-84坐标系 srs = osr.SpatialReference() srs.ImportFromEPSG(4326) # WGS-84的EPSG代码 dst_ds.SetProjection(srs.ExportToWkt()) dst_ds.FlushCache() return True

4. 处理图像旋转问题

大疆无人机在拍摄时可能会旋转角度，这会导致照片方向与实际地理方向不一致。我们需要从元数据中获取偏航角（Yaw）并进行校正。

关键步骤：

1. 从drone-dji:FlightYawDegree字段获取旋转角度
2. 计算图像中心点坐标
3. 应用仿射变换进行旋转

旋转校正的核心算法：

from affine import Affine def get_yaw_angle(image_path): """从图像元数据中提取偏航角""" with open(image_path, 'rb') as f: content = f.read().decode('latin-1') match = re.search(r'drone-dji:FlightYawDegree="([^"]+)"', content) if match: return float(match.group(1)) return 0.0 def apply_rotation(input_path, output_path): yaw = get_yaw_angle(input_path) if yaw == 0: return False # 无需旋转 src_ds = gdal.Open(input_path) gt = src_ds.GetGeoTransform() # 计算中心点坐标 x_size, y_size = src_ds.RasterXSize, src_ds.RasterYSize center_x = gt[0] + x_size/2 * gt[1] center_y = gt[3] + y_size/2 * gt[5] # 创建旋转后的仿射变换 affine_src = Affine.from_gdal(*gt) affine_dst = affine_src * affine_src.rotation(yaw, (center_x, center_y)) # 保存结果 driver = gdal.GetDriverByName('GTiff') dst_ds = driver.CreateCopy(output_path, src_ds) dst_ds.SetGeoTransform(affine_dst.to_gdal()) dst_ds.FlushCache() return True

5. 完整工作流与实战示例

现在我们将上述步骤整合成一个完整的处理流程，并提供一个实际应用案例。

典型工作流程：

1. 输入原始无人机照片
2. 提取GPS位置和图像尺寸信息
3. 计算仿射变换参数
4. 应用WGS-84坐标系
5. 校正图像旋转（如果需要）
6. 输出带有地理参考的TIFF文件

完整代码示例：

import os from osgeo import gdal, osr import re from affine import Affine class DJI_GeoReference: def __init__(self, resolution=0.00001): self.resolution = resolution # 默认分辨率（度/像素） def process_image(self, input_path, output_folder): """处理单张无人机照片""" if not os.path.exists(input_path): raise FileNotFoundError(f"输入文件不存在: {input_path}") os.makedirs(output_folder, exist_ok=True) base_name = os.path.splitext(os.path.basename(input_path))[0] temp_path = os.path.join(output_folder, f"{base_name}_temp.tif") output_path = os.path.join(output_folder, f"{base_name}_georef.tif") # 第一步：添加坐标系 self._assign_coordinate_system(input_path, temp_path) # 第二步：校正旋转 if self._apply_rotation(temp_path, output_path): os.remove(temp_path) # 删除临时文件 else: os.rename(temp_path, output_path) return output_path def _extract_metadata(self, image_path): """提取元数据的内部方法""" # ...（同前面的extract_metadata函数） def _create_geotransform(self, metadata): """创建仿射变换的内部方法""" # ...（同前面的create_geotransform函数） def _assign_coordinate_system(self, input_path, output_path): """添加坐标系的内部方法""" # ...（同前面的assign_coordinate_system函数） def _get_yaw_angle(self, image_path): """获取偏航角的内部方法""" # ...（同前面的get_yaw_angle函数） def _apply_rotation(self, input_path, output_path): """应用旋转的内部方法""" # ...（同前面的apply_rotation函数） # 使用示例 if __name__ == "__main__": processor = DJI_GeoReference(resolution=0.00001) input_image = "DJI_001.JPG" # 替换为你的照片路径 output_dir = "output_georef" result = processor.process_image(input_image, output_dir) print(f"处理完成，结果保存在: {result}")

注意：实际分辨率取决于无人机飞行高度和相机参数。对于大疆御3E，在100米高度拍摄时，典型分辨率约为3厘米/像素，对应的经纬度分辨率大约为0.00001度/像素。

6. 验证与调试技巧

处理后的影像需要验证其地理定位的准确性。以下是几种实用的验证方法：

• QGIS可视化检查：将处理后的影像加载到QGIS中，与已知的基准地图（如OpenStreetMap）叠加，检查位置是否匹配
• 控制点验证：在地面布设已知坐标的控制点，检查影像上对应点的坐标
• 元数据交叉验证：比较处理后的影像元数据与原始EXIF中的GPS信息

常见问题及解决方案：

调试时可以添加详细的日志输出：

import logging logging.basicConfig(level=logging.INFO) logger = logging.getLogger(__name__) def debug_geo_transform(dataset): """打印详细的坐标变换信息""" gt = dataset.GetGeoTransform() logger.info(f"仿射变换参数: {gt}") logger.info(f"左上角坐标: ({gt[0]}, {gt[3]})") logger.info(f"像素大小: X={gt[1]}度, Y={gt[5]}度") logger.info(f"旋转参数: X={gt[2]}, Y={gt[4]}")

7. 性能优化与批量处理

当需要处理大量无人机照片时，效率变得尤为重要。以下是几种优化策略：

并行处理技术：

from concurrent.futures import ThreadPoolExecutor def batch_process(image_list, output_dir, workers=4): """批量处理无人机照片""" processor = DJI_GeoReference() with ThreadPoolExecutor(max_workers=workers) as executor: futures = [] for img_path in image_list: future = executor.submit( processor.process_image, img_path, output_dir ) futures.append(future) results = [f.result() for f in futures] return results

内存优化技巧：

• 使用分块处理大图像
• 及时释放GDAL数据集对象
• 适当降低输出图像的质量等级以减小文件大小

def memory_efficient_process(input_path, output_path): """内存优化的处理方式""" # 创建输出文件时设置压缩选项 creation_options = [ 'COMPRESS=LZW', 'PREDICTOR=2', 'TILED=YES' ] src_ds = gdal.Open(input_path) driver = gdal.GetDriverByName('GTiff') dst_ds = driver.CreateCopy( output_path, src_ds, options=creation_options ) # ...（其余处理逻辑） # 显式释放资源 src_ds = None dst_ds = None

在实际项目中，处理100张大疆御3E照片（每张约20MB）的典型性能指标：