首页 > 编程语言 >【Python&RS】基于GDAL给无人机图片定义坐标系

【Python&RS】基于GDAL给无人机图片定义坐标系

时间：2023-07-13 10:45:17浏览次数：40  
标签：angle RS Python sum meter height GDAL ds math

​

        前段时间有过一个想法，就是如果可以给无人机拍摄的图片定义坐标系，再使用GADL库里的镶嵌拼接函数，是不是就可以实现快速拼接影像。虽然结果不是正射影像，但效率比无人机厂家的软件提高了很多很多，主要还是看用途。

        有了这个想法后就要行动起来，定义一个坐标系，一般需要两个点。一是需要投影坐标系的参数，二就是仿射地理变换参数。投影坐标系非常容易得到，投影坐标系的参数GDAL有内置的一部分，当然也可以读取已有数据的坐标系。但仿射地理变换参数就比较麻烦了，我印象中定义仿射地理变换参数时只能通过图片左上角的点就是设置，不像ENVI可以任意选点。ENVI方便很多，有兴趣可以查看【RS】基于ENVI给图片/影像定义坐标系。所以我们只要能解决仿射地理变换参数就可以解决图片定义坐标系的问题。

​​

 

大致思路是：通过之前计算图片某点的投影坐标的方法，得到图片左上角的投影坐标和偏转角度，以此建立仿射地理变换参数，再通过这些创建新的数据资源将原来的数据写入。

 


一、获取图片中心点经纬度

        这一步之前的博文已经说过无数次了，这里就不解释。如果你只是娱乐一下，并不需要定义真实的坐标系，那么可以直接跳到最后一步，仿射地理变换参数直接设置0，1即可。

def Get_LatLon(path_image):
    """
    :param path_image: 输入图片路径
    :return: 返回图片中心点经纬度
    """
    # 获取图片的经纬度信息
    f = open(path_image, 'rb')
    contents = exifread.process_file(f)
    longitude = contents["GPS GPSLongitude"].values
    longitude_f = longitude[0].num/longitude[0].den + (longitude[1].num/longitude[1].den/60) + (longitude[2].num/longitude[2].den/3600)
    latitude = contents["GPS GPSLatitude"].values
    latitude_f = latitude[0].num/latitude[0].den + (latitude[1].num/latitude[1].den/60) + (latitude[2].num/latitude[2].den/3600)
    # print("经度:", longitude_f)  # contents['GPS GPSLatitudeRef']
    # print("纬度:", latitude_f)  # contents['GPS GPSLongitudeRef']
    f.close()
    return latitude_f, longitude_f

二、转换投影

        无人机拍摄的影像一般都是GPS经纬度坐标（WGS84），我这里要将他定义为UTM/WGS84 51N，所以需要将它的地理坐标系先转换成投影坐标系。

def LonLat_Meter(lat, lon):
    """
    :param lat: 图片中心点纬度
    :param lon: 图片中心点经度
    :return: 返回UTM投影坐标
    """
    source = osr.SpatialReference()
    # 初始化osr.SpatialReference对象以形成一个合法的坐标系统
    source.ImportFromEPSG(4326)
    # 向对象中写入Source_EPSG坐标系统
    target = osr.SpatialReference()
    target.ImportFromEPSG(32651)
    # 这里是用内置的EPSG对应的坐标系作为转换参数
    transform = osr.CoordinateTransformation(source, target)
    point = ogr.CreateGeometryFromWkt("POINT (%s %s)" % (lat, lon))
    # 报错的话，将经纬度翻过来
    point.Transform(transform)
    # print(point.GetX(), point.GetY())
    return point.GetX(), point.GetY()

三、获取图片偏转角

        无人机航拍时不可能一直与正北方向一致，所以还需要获取其成像时的偏转角。这一步是为了计算图片左上角的仿射地理变换参数。我这里使用的是大疆的无人机，其他厂商的无人机成像参数需要自己查阅。

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 Target_Meter(Image_Yaw_angle, width, height, meter_X, meter_Y, rat):
    """
    :param Image_Yaw_angle: 图片与正北的偏转角
    :param width: 图片的宽度的一半
    :param height: 图片的高度的一半
    :param meter_X: 图片中心点的投影坐标
    :param meter_Y: 图片中心点的投影坐标
    :param rat: 图片分辨率
    :return: 返回左上角的投影坐标
    """
    if Image_Yaw_angle < 0:
        Image_Yaw_angle = 360 + Image_Yaw_angle
    yaw_angle = 360 - math.degrees(math.atan(width / height))
    # 左上角与图片正北的夹角
    sum_angle = Image_Yaw_angle + yaw_angle
    # 图片偏转角+左上角偏转角
    if sum_angle >= 360:
        sum_angle = sum_angle - 360
    if sum_angle == 0:
        meter_x = meter_X
        meter_y = meter_Y - height*rat
    elif sum_angle == 90:
        meter_x = meter_X + (-width)*rat
        meter_y = meter_Y
    elif sum_angle == 180:
        meter_x = meter_X
        meter_y = meter_Y + (-height)*rat
    elif sum_angle == 270:
        meter_x = meter_X - width*rat
        meter_y = meter_Y
    elif sum_angle == 360:
        meter_x = meter_X
        meter_y = meter_Y - height*rat
    elif 0 < sum_angle < 90:
        meter_x = meter_X + math.sin(math.radians(sum_angle))*math.sqrt(math.pow(-width, 2)+math.pow(-height, 2))*rat
        meter_y = meter_Y + math.cos(math.radians(sum_angle))*math.sqrt(math.pow(-width, 2)+math.pow(-height, 2))*rat
    elif 90 < sum_angle < 180:
        meter_x = meter_X + math.sin(math.radians(180-sum_angle)) * math.sqrt(math.pow(-width, 2) +
                                                                              math.pow(-height, 2))*rat
        meter_y = meter_Y - math.cos(math.radians(180-sum_angle)) * math.sqrt(math.pow(-width, 2) +
                                                                              math.pow(-height, 2))*rat
    elif 180 < sum_angle < 270:
        meter_x = meter_X - math.sin(math.radians(sum_angle-180)) * math.sqrt(math.pow(-width, 2) +
                                                                              math.pow(-height, 2))*rat
        meter_y = meter_Y - math.cos(math.radians(sum_angle-180)) * math.sqrt(math.pow(-width, 2) +
                                                                              math.pow(-height, 2))*rat
    elif 270 < sum_angle < 360:
        meter_x = meter_X - math.sin(math.radians(360-sum_angle)) * math.sqrt(math.pow(-width, 2) +
                                                                              math.pow(-height, 2))*rat
        meter_y = meter_Y + math.cos(math.radians(360-sum_angle)) * math.sqrt(math.pow(-width, 2) +
                                                                              math.pow(-height, 2))*rat
    return meter_x, meter_y

五、定义坐标系

        前面的工作完成后就得到了仿射地理变换参数，就可以定义坐标系了。如果大家只是想拼接手机拍摄的图片，或者娱乐学习，就可以不那么严谨。不需要求出真实的投影坐标，仿射地理变换参数直接设置为(1, 0, 0, 1, 0, 0)即可。关于仿射地理变换参数的含义，大家可以查看这篇博文【Python&RS】遥感影像的像素坐标转地理坐标（仿射变换）。

def Define_Proj(input_path, ds_geo, out_path):
    """
    :param input_path: 需要定义投影的图片路径
    :param ds_geo: 仿射地理变换参数，元组的形式
    :param out_path: 输出路径
    :return:
    """
    ds = gdal.Open(input_path)  # 打开数据集dataset，需要定义投影的图片
    ds_width = ds.RasterXSize  # 获取数据宽度
    ds_height = ds.RasterYSize  # 获取数据高度
    ds_bands = ds.RasterCount  # 获取波段数
    driver = gdal.GetDriverByName('GTiff')  # 载入数据驱动，用于存储内存中的数组
    ds_result = driver.Create(out_path, ds_width, ds_height, bands=ds_bands, eType=gdal.GDT_Float64)
    ds_result.SetGeoTransform(ds_geo)  # 导入仿射地理变换参数
    target = osr.SpatialReference()
    target.ImportFromEPSG(32651)
    # 输出投影坐标系为WGS84 UTM 51N
    ds_result.SetProjection(str(target))  # 导入投影信息
    for i in range(1, ds_bands+1):
        band_data = ds.GetRasterBand(i).ReadAsArray(0, 0, ds_width, ds_height).astype(np.cfloat)
        ds_result.GetRasterBand(i).WriteArray(band_data)
    del ds_result

六、完整代码

        上面的代码都是分步的，直接复制无法运行，因为有一部分我是放在了主函数里面。运行代码要注意两点，一是无人机图片是否有中心点的经纬度，二是是否为大疆的无人机，如果不是第三步的代码就需要自己修改。

# -*- coding: utf-8 -*-
"""
@Time ： 2023/5/29 15:28
@Auth ： RS迷途小书童
@File ：Define Coordinate System For Image.py
@IDE ：PyCharm
@Purpose ：给图片定义坐标系
"""
import exifread
from osgeo import gdal, ogr, osr
import math
import numpy as np


def Get_LatLon(path_image):
    """
    :param path_image: 输入图片路径
    :return: 返回图片中心点经纬度
    """
    # 获取图片的经纬度信息
    f = open(path_image, 'rb')
    contents = exifread.process_file(f)
    longitude = contents["GPS GPSLongitude"].values
    longitude_f = longitude[0].num/longitude[0].den + (longitude[1].num/longitude[1].den/60) + (longitude[2].num/longitude[2].den/3600)
    latitude = contents["GPS GPSLatitude"].values
    latitude_f = latitude[0].num/latitude[0].den + (latitude[1].num/latitude[1].den/60) + (latitude[2].num/latitude[2].den/3600)
    # print("经度:", longitude_f)  # contents['GPS GPSLatitudeRef']
    # print("纬度:", latitude_f)  # contents['GPS GPSLongitudeRef']
    f.close()
    return latitude_f, longitude_f


def LonLat_Meter(lat, lon):
    """
    :param lat: 图片中心点纬度
    :param lon: 图片中心点经度
    :return: 返回UTM投影坐标
    """
    source = osr.SpatialReference()
    # 初始化osr.SpatialReference对象以形成一个合法的坐标系统
    source.ImportFromEPSG(4326)
    # 向对象中写入Source_EPSG坐标系统
    target = osr.SpatialReference()
    target.ImportFromEPSG(32651)
    # 这里是用内置的EPSG对应的坐标系作为转换参数
    transform = osr.CoordinateTransformation(source, target)
    point = ogr.CreateGeometryFromWkt("POINT (%s %s)" % (lat, lon))
    # 报错的话，将经纬度翻过来
    point.Transform(transform)
    # print(point.GetX(), point.GetY())
    return point.GetX(), point.GetY()


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 Target_Meter(Image_Yaw_angle, width, height, meter_X, meter_Y, rat):
    """
    :param Image_Yaw_angle: 图片与正北的偏转角
    :param width: 图片的宽度的一半
    :param height: 图片的高度的一半
    :param meter_X: 图片中心点的投影坐标
    :param meter_Y: 图片中心点的投影坐标
    :param rat: 图片分辨率
    :return: 返回左上角的投影坐标
    """
    if Image_Yaw_angle < 0:
        Image_Yaw_angle = 360 + Image_Yaw_angle
    yaw_angle = 360 - math.degrees(math.atan(width / height))
    # 左上角与图片正北的夹角
    sum_angle = Image_Yaw_angle + yaw_angle
    # 图片偏转角+左上角偏转角
    if sum_angle >= 360:
        sum_angle = sum_angle - 360
    if sum_angle == 0:
        meter_x = meter_X
        meter_y = meter_Y - height*rat
    elif sum_angle == 90:
        meter_x = meter_X + (-width)*rat
        meter_y = meter_Y
    elif sum_angle == 180:
        meter_x = meter_X
        meter_y = meter_Y + (-height)*rat
    elif sum_angle == 270:
        meter_x = meter_X - width*rat
        meter_y = meter_Y
    elif sum_angle == 360:
        meter_x = meter_X
        meter_y = meter_Y - height*rat
    elif 0 < sum_angle < 90:
        meter_x = meter_X + math.sin(math.radians(sum_angle))*math.sqrt(math.pow(-width, 2)+math.pow(-height, 2))*rat
        meter_y = meter_Y + math.cos(math.radians(sum_angle))*math.sqrt(math.pow(-width, 2)+math.pow(-height, 2))*rat
    elif 90 < sum_angle < 180:
        meter_x = meter_X + math.sin(math.radians(180-sum_angle)) * math.sqrt(math.pow(-width, 2) +
                                                                              math.pow(-height, 2))*rat
        meter_y = meter_Y - math.cos(math.radians(180-sum_angle)) * math.sqrt(math.pow(-width, 2) +
                                                                              math.pow(-height, 2))*rat
    elif 180 < sum_angle < 270:
        meter_x = meter_X - math.sin(math.radians(sum_angle-180)) * math.sqrt(math.pow(-width, 2) +
                                                                              math.pow(-height, 2))*rat
        meter_y = meter_Y - math.cos(math.radians(sum_angle-180)) * math.sqrt(math.pow(-width, 2) +
                                                                              math.pow(-height, 2))*rat
    elif 270 < sum_angle < 360:
        meter_x = meter_X - math.sin(math.radians(360-sum_angle)) * math.sqrt(math.pow(-width, 2) +
                                                                              math.pow(-height, 2))*rat
        meter_y = meter_Y + math.cos(math.radians(360-sum_angle)) * math.sqrt(math.pow(-width, 2) +
                                                                              math.pow(-height, 2))*rat
    return meter_x, meter_y


def Define_Proj(input_path, ds_geo, out_path):
    """
    :param input_path: 需要定义投影的图片路径
    :param ds_geo: 仿射地理变换参数，元组的形式
    :param out_path: 输出路径
    :return:
    """
    ds = gdal.Open(input_path)  # 打开数据集dataset，需要定义投影的图片
    ds_width = ds.RasterXSize  # 获取数据宽度
    ds_height = ds.RasterYSize  # 获取数据高度
    ds_bands = ds.RasterCount  # 获取波段数
    driver = gdal.GetDriverByName('GTiff')  # 载入数据驱动，用于存储内存中的数组
    ds_result = driver.Create(out_path, ds_width, ds_height, bands=ds_bands, eType=gdal.GDT_Float64)
    ds_result.SetGeoTransform(ds_geo)  # 导入仿射地理变换参数
    target = osr.SpatialReference()
    target.ImportFromEPSG(32651)
    # 输出投影坐标系为WGS84 UTM 51N
    ds_result.SetProjection(str(target))  # 导入投影信息
    for i in range(1, ds_bands+1):
        band_data = ds.GetRasterBand(i).ReadAsArray(0, 0, ds_width, ds_height).astype(np.cfloat)
        ds_result.GetRasterBand(i).WriteArray(band_data)
    del ds_result


if __name__ == '__main__':
    path_image1 = 'G:/A_01_DSC00005.JPG'
    out_image = 'G:/try/3.tif'
    latitude1, longitude1 = Get_LatLon(path_image1)
    # 获取图片中心点的纬度、经度
    X, Y = LonLat_Meter(latitude1, longitude1)
    x, y = Target_Meter(135, 3000, 2000, X, Y, 0.04)
    geo = tuple((x, 0.04, 0, y, 0, 0.04))
    print(geo)
    Define_Proj(path_image1, geo, out_image)

 

        我这里只是实验一下，所以仿射地理变换参数里的旋转角度没有加进去，位置有一点偏移。总体代码没有问题，但这个逻辑应该还需要改进。大家有好的想法也可以和我一起探讨，后续应该会针对这个问题进行优化。

        前端时间去出外业了，所以博文没有更新。后面应该会持续更新的，当然博主也会休息休息，偷偷懒，要一起加油鸭。

​

标签：angle,RS,Python,sum,meter,height,GDAL,ds,math
From： https://www.cnblogs.com/RSran/p/17549764.html

相关文章

  • 【Python】数据可视化利器PyCharts在测试工作中的应用
    PyCharts简介PyCharts是一个基于Python的数据可视化库，它支持多种图表类型，如折线图、柱状图、饼图等。PyCharts提供了简洁的API，使得用户能够轻松地创建各种图表，同时支持个性化的配置，以满足不同需求。PyCharts的底层依赖于ECharts，这使得它在功能和性能上都具有很高的优势。......
  • ubuntu20使用iptables-persistent libfakeroot libxtables-dev netfilter-persistent
    实施防火墙是保护服务器安全的重要一步。其中很大一部分是决定将对您的网络实施流量限制的单个规则和策略。像iptables这样的防火墙还允许您对应用规则的结构框架有发言权。在本指南中,您将学习如何构建防火墙,作为更复杂规则集的基础。该防火墙将主要关注提供合理的默认值和建立......
  • PYTHON随笔-logging
    PYTHON随笔-loggingimportloggingfromlogging.handlersimportRotatingFileHandlergLogFile='/home/mcs/log/dbm_py.log'LOG_FORMAT="%(asctime)s[%(levelname)-5s][%(filename)s:%(lineno)s]-%(message)s"#日志输出的格式#-8表示占位符，让输出左对齐，输......
  • python-pymysql-类对象映射为sql语句
    查询语句importpymysqlclassUserQuery:def__init__(self,name=None,age=None,email=None):self.name=nameself.age=ageself.email=emaildefselect_data(table,condition):#连接到数据库connection=pymysql.connec......
  • dircolors
    dircolors置ls命令在显示目录或文件时所用的色彩补充说明dircolors命令设置ls命令在显示目录或文件时所用的色彩。dircolors可根据[色彩配置文件]来设置LS_COLORS环境变量或是显示设置LS_COLORS环境变量的命令。语法dircolors(选项)(参数)选项-b或--sh或--bourne-shell：显......
  • Python用Keras神经网络序列模型回归拟合预测、准确度检查和结果可视化|附代码数据
    原文链接：http://tecdat.cn/?p=23573最近我们被客户要求撰写关于Keras神经网络序列模型的研究报告，包括一些图形和统计输出。我们可以很容易地用Keras序列模型拟合回归数据并预测测试数据。  在这篇文章中，我们将简要地学习如何用Python中的Keras神经网络API拟合回归数据。我们将......
  • 双服务台串联排队系统——Python仿真
    串联排队系统是一种常见的排队系统结构，由多个单一排队系统按照特定的顺序连接而成。在串联排队系统中，一个排队系统的输出将成为下一个排队系统的输入，从而形成连续的流程。这种系统结构可以用于模拟和优化许多现实世界中的流程，如生产线、物流运输等。一、双服务台串联排队系统模......
  • 如何使用Python制作交互式股票K线图？
    如何使用Python制作交互式股票K线图？如何使用Python制作交互式股票K线图？-知乎(zhihu.com)州的先生  在之前的文章中，我们介绍了使用PyQtGraph在PyQt5中绘制股票K线图：PythonGUI教程（十三）：在GUI中使用pyqtgraph绘图库​zmister.com/archives/187.html以及使......
  • Python3.6下scrapy框架的安装
    命令安装，提示  FailedbuildingwheelforTwistedMicrosoftVisualC++14.0isrequired...  总结pipinstallwheel 下载Twisted包安装下载Scrapy包安装下载地址：http://www.lfd.uci.edu/~gohlke/pythonlibs/详细解决方案1首先考虑使用最简单的方法安装pipinstallsc......
  • Python-[]列表.py
     19printlist;            #输出完整列表 20printlist[0]  #输出列表第一个元素 21printlist[1:3]#输出列表下标1～3之间的元素（和字符串一样，含头不含尾） 22printlist[2:] #输出下标2以后所有的元素（包含下标2的元素） 23printtinylist*2     ......

赞助商

阅读排行

网站主页关于我们联系我们网站地图

本网站内容转载自其他媒体，侵权联系[admin##ips99.com]。

Copyright © 2020-2023 IPS99 版权所有 IPS99