文章详情

短信预约-IT技能 免费直播动态提醒

请输入下面的图形验证码

提交验证

短信预约提醒成功

【Python&GIS】无人机影像的像素坐标计算图片某点的地理/投影坐标

2023-09-03 09:13

关注

        又是掉头发的一天,今天的任务是通过图片中心点的地理坐标以及图片中某点的像素坐标(即这个点位于图片中的x,y坐标)计算该点的地理/投影坐标。经过一整天的搜索,发现网上并没有这方面的教程。然后就是想啊想,头发一抓一大把,终于在网上零零散散的教程以及不断摸索中解决了这个问题。

        大致思路就是,先获取图片相对真北方向的偏转角以及该点和图片中心的连线与图片的正北方向夹角;然后将图片中心点的地理坐标转换为投影坐标(如果这一步没有中心点的地理坐标,那么你就不用继续往下看了);最后就是通过图片分辨率计算点到中心的实际距离,再通过夹角和中心点的投影坐标加加减减即可。话虽这么说,但实施起来真心不简单啊!!!

一、准备工作

1.获取图片中心点的经纬度坐标/投影坐标(必须)

        如果只有一张图片的话,你完全可以右键图片,查看其属性,里面就有经纬度坐标。同时也可以使用Python实现,之前分享过【Python&GIS】判断图片中心点/经纬度点是否在某个面内,所以这里就不解释了,直接上代码。

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)    f.close()    return longitude_f, latitude_f

2.获取图片与真北方向的偏转角(非必须)

        如果你的图片并不是无人机拍的,或者拍摄时图片与真北方向并无夹角,那么就不需要这一步。这一步主要就是矫正相机拍摄时的偏转,每种无人机的参数名可能不一样,所以需要自行修改。

        下面的代码也可以获取图片中心点的经纬度,但我之前以及用过第一步的代码了,所以就继续使用第一步的了,你们可以自己简化。

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])

3.获取图片目标点与中心点的连线和图片正北方向的夹角(必须)

        通俗来说,就是该点与图片正上方的夹角,同样这步也是用来矫正偏转的,获取该角度后,即可得到该点在投影坐标系中与真北方向的夹角。

def Yaw_angle_calculation(x, y, x_Center, y_Center):    """    :param x: 输入图片目标点的栅格坐标x    :param y: 输入图片目标点的栅格坐标y    :param x_Center: 输入图片的中心点x坐标    :param y_Center: 输入图片的中心点y坐标    :return: 目标点相对于中心点的偏转角度    """    # 获取目标的偏转角    if x == x_Center and y == y_Center:        Yaw_angle = 0    elif x == x_Center and y != y_Center:        if y > y_Center:            Yaw_angle = 180        elif y < y_Center:            Yaw_angle = 0    elif x != x_Center and y == y_Center:        if x >x_Center:            Yaw_angle = 90        elif x < x_Center:            Yaw_angle = 270    else:        if x > x_Center:            if y < y_Center:                Yaw_angle = math.degrees(math.atan((x - x_Center) / (y_Center - y)))            else:                Yaw_angle = 180 - math.degrees(math.atan((x - x_Center) / (y - y_Center)))        else:            if y < y_Center:                Yaw_angle = 360 - math.degrees(math.atan((x_Center - x) / (y_Center - y)))            else:                Yaw_angle = 180 + math.degrees(math.atan((x_Center - x) / (y - y_Center)))    # print("Yaw_angle",Yaw_angle)    return Yaw_angle

4.将图片中心点的地理坐标转为投影坐标(必须)

        一般来说图片的中心点都是GPS坐标,即WGS84地理坐标系的经纬度。但当我们计算距离时,需要使用到投影坐标系,所以这里需要将其转换一下。我这里是WGS84地理坐标系转为UTM   51N投影坐标系,你们视情况修改。

def LonLat_Meter(lat, lon):    """    :param lat: 输入WGS84经度    :param lon: 输入WGS84纬度    :return: 返回WGS84/UTM 51N的x,y    """    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()

二、转换函数

1.前期工作完成后,就可以实现转换了(必须!)

        这里输入参数为图片偏转角、目标点的像素坐标(栅格坐标)、中心点的像素坐标、栅格分辨率(空间分辨率)、中心点的投影坐标。

def Target_Meter(Image_Yaw_angle, x, y, x_Center, y_Center, ratio,  meter_X, meter_Y):    """    :param Image_Yaw_angle: 输入图片的偏航角    :param x: 输入图片目标点的栅格坐标x    :param y: 输入图片目标点的栅格坐标y    :param x_Center: 输入图片中心点的栅格坐标x    :param y_Center: 输入图片中心点的栅格坐标y    :param meter_X: 输入图片中心点的投影坐标x(UTM/WGS84 51N)    :param meter_Y: 输入图片中心点的投影坐标y(UTM/WGS84 51N)    :return: 目标点的投影坐标x,y    """    if Image_Yaw_angle < 0:        Image_Yaw_angle = 360 +Image_Yaw_angle    Yaw_angle = Yaw_angle_calculation(x, y, x_Center, y_Center)    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 - (y_Center-y)*ratio    elif sum_angle == 90:        meter_x = meter_X + (x-x_Center)*ratio        meter_y = meter_Y    elif sum_angle == 180:        meter_x = meter_X        meter_y = meter_Y + (y-y_Center)*ratio    elif sum_angle == 270:        meter_x = meter_X - (x_Center-x)*ratio        meter_y = meter_Y    elif sum_angle == 360:        meter_x = meter_X        meter_y = meter_Y - (x_Center-y)*ratio    elif 0 < sum_angle < 90:        meter_x = meter_X + math.sin(math.radians(sum_angle))*math.sqrt(math.pow(x-x_Center, 2)+math.pow(y-y_Center, 2))*ratio        meter_y = meter_Y + math.cos(math.radians(sum_angle))*math.sqrt(math.pow(x-x_Center, 2)+math.pow(y-y_Center, 2))*ratio    elif 90 < sum_angle < 180:        meter_x = meter_X + math.sin(math.radians(180-sum_angle)) * math.sqrt(math.pow(x - x_Center, 2) + math.pow(y -y_Center, 2))*ratio        meter_y = meter_Y - math.cos(math.radians(180-sum_angle)) * math.sqrt(math.pow(x - x_Center, 2) + math.pow(y -y_Center, 2))*ratio    elif 180 < sum_angle < 270:        meter_x = meter_X - math.sin(math.radians(sum_angle-180)) * math.sqrt(math.pow(x - x_Center, 2) + math.pow(y - y_Center, 2))*ratio        meter_y = meter_Y - math.cos(math.radians(sum_angle-180)) * math.sqrt(math.pow(x - x_Center, 2) + math.pow(y - y_Center, 2))*ratio    elif 270 < sum_angle < 360:        meter_x = meter_X - math.sin(math.radians(360-sum_angle)) * math.sqrt(math.pow(x - x_Center, 2) + math.pow(y - y_Center, 2))*ratio        meter_y = meter_Y + math.cos(math.radians(360-sum_angle)) * math.sqrt(math.pow(x - x_Center, 2) + math.pow(y - y_Center, 2))*ratio    return meter_x, meter_y

2.将目标点的投影坐标转为地理坐标(非必须)

        我这里因为项目最后需要的是经纬度坐标,所以在计算得到目标点的投影坐标后,还需要转换一下。你们视情况而定,如果你们只要投影坐标,那么第5步输出的就已经是投影坐标啦。   

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

 三、完整代码

# -*- coding: utf-8 -*-"""@Time : 2023/1/5 17:45@Auth : RS迷途小书童@File :Picture coordinate system to geographic coordinate.py@IDE :PyCharm@Purpose :图片某点的像素坐标转为地理坐标系"""import mathimport exifreadfrom osgeo import ogr, osrdef 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)    f.close()    return longitude_f, latitude_fdef LonLat_Meter(lat, lon):    """    :param lat: 输入WGS84经度    :param lon: 输入WGS84纬度    :return: 返回WGS84/UTM 51N的x,y    """    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 Meter_LonLat(x, y):    """    :param x: 输入WGS84/UTM 51N的x    :param y: 输入WGS84/UTM 51N的y    :return: 返回WGS84的经纬度    """    source = osr.SpatialReference()    # 初始化osr.SpatialReference对象以形成一个合法的坐标系统    source.ImportFromEPSG(32651)    # 向对象中写入Source_EPSG坐标系统    target = osr.SpatialReference()    target.ImportFromEPSG(4326)    # 这里是用内置的EPSG对应的坐标系作为转换参数    transform = osr.CoordinateTransformation(source, target)    point = ogr.CreateGeometryFromWkt("POINT (%s %s)" % (x, y))    # 报错的话,将经纬度翻过来    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 Yaw_angle_calculation(x, y, x_Center, y_Center):    """    :param x: 输入图片目标点的栅格坐标x    :param y: 输入图片目标点的栅格坐标y    :return: 目标点相对于中心点的偏转角度    """    # 获取目标的偏转角    if x == x_Center and y == y_Center:        Yaw_angle = 0    elif x == x_Center and y != y_Center:        if y > y_Center:            Yaw_angle = 180        elif y < 1500:            Yaw_angle = 0    elif x != x_Center and y == y_Center:        if x >x_Center:            Yaw_angle = 90        elif x < x_Center:            Yaw_angle = 270    else:        if x > x_Center:            if y < y_Center:                Yaw_angle = math.degrees(math.atan((x - x_Center) / (y_Center - y)))            else:                Yaw_angle = 180 - math.degrees(math.atan((x - x_Center) / (y - y_Center)))        else:            if y < y_Center:                Yaw_angle = 360 - math.degrees(math.atan((x_Center - x) / (y_Center - y)))            else:                Yaw_angle = 180 + math.degrees(math.atan((x_Center - x) / (y - y_Center)))    # print("Yaw_angle",Yaw_angle)    return Yaw_angledef Target_Meter(Image_Yaw_angle, x, y, x_Center, y_Center, ratio,  meter_X, meter_Y):    """    :param Image_Yaw_angle: 输入图片的偏航角    :param x: 输入图片目标点的栅格坐标x    :param y: 输入图片目标点的栅格坐标y    :param x_Center: 输入图片中心点的栅格坐标x    :param y_Center: 输入图片中心点的栅格坐标y    :param meter_X: 输入图片中心点的投影坐标x(UTM/WGS84 51N)    :param meter_Y: 输入图片中心点的投影坐标y(UTM/WGS84 51N)    :return: 目标点的投影坐标x,y    """    if Image_Yaw_angle < 0:        Image_Yaw_angle = 360 +Image_Yaw_angle    Yaw_angle = Yaw_angle_calculation(x, y, x_Center, y_Center)    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 - (y_Center-y)*ratio    elif sum_angle == 90:        meter_x = meter_X + (x-x_Center)*ratio        meter_y = meter_Y    elif sum_angle == 180:        meter_x = meter_X        meter_y = meter_Y + (y-y_Center)*ratio    elif sum_angle == 270:        meter_x = meter_X - (x_Center-x)*ratio        meter_y = meter_Y    elif sum_angle == 360:        meter_x = meter_X        meter_y = meter_Y - (x_Center-y)*ratio    elif 0 < sum_angle < 90:        meter_x = meter_X + math.sin(math.radians(sum_angle))*math.sqrt(math.pow(x-x_Center, 2)+math.pow(y-y_Center, 2))*ratio        meter_y = meter_Y + math.cos(math.radians(sum_angle))*math.sqrt(math.pow(x-x_Center, 2)+math.pow(y-y_Center, 2))*ratio    elif 90 < sum_angle < 180:        meter_x = meter_X + math.sin(math.radians(180-sum_angle)) * math.sqrt(math.pow(x - x_Center, 2) + math.pow(y -y_Center, 2))*ratio        meter_y = meter_Y - math.cos(math.radians(180-sum_angle)) * math.sqrt(math.pow(x - x_Center, 2) + math.pow(y -y_Center, 2))*ratio    elif 180 < sum_angle < 270:        meter_x = meter_X - math.sin(math.radians(sum_angle-180)) * math.sqrt(math.pow(x - x_Center, 2) + math.pow(y - y_Center, 2))*ratio        meter_y = meter_Y - math.cos(math.radians(sum_angle-180)) * math.sqrt(math.pow(x - x_Center, 2) + math.pow(y - y_Center, 2))*ratio    elif 270 < sum_angle < 360:        meter_x = meter_X - math.sin(math.radians(360-sum_angle)) * math.sqrt(math.pow(x - x_Center, 2) + math.pow(y - y_Center, 2))*ratio        meter_y = meter_Y + math.cos(math.radians(360-sum_angle)) * math.sqrt(math.pow(x - x_Center, 2) + math.pow(y - y_Center, 2))*ratio    return meter_x, meter_yif __name__ == "__main__":    path = "G:/DJI_26_W.jpeg"    x = 1600    y = 2500    x_Center = 2000    y_Center = 1500    ratio = 0.046    longitude_f, latitude_f = Get_LatLon(path)    meter_X, meter_Y = LonLat_Meter(latitude_f, longitude_f)    Image_Yaw_angle = Get_Image_Yaw_angle(path)    # 获取图片的偏转角    meter_x, meter_y = Target_Meter(Image_Yaw_angle, x, y, x_Center, y_Center, ratio,  meter_X, meter_Y)    lat_x, lon_y = Meter_LonLat(meter_x, meter_y)    print(longitude_f, latitude_f)    print(meter_X, meter_Y)    print(Image_Yaw_angle)    print(meter_x, meter_y)    print(lat_x, lon_y)

        大家可以自己添加for、while循环实现多张图片或图片中多个目标物的坐标转换。博主这里是为了目标识别后的目标物真实坐标的计算,所以将循环部分内置到了其他程序。

        本文章主要是分享个人在学习Python过程中写过的一些代码。有些部分参考了前人以及官网的教程,如有侵权请联系作者删除,大家有问题可以随时留言交流,博主会及时回复。

来源地址:https://blog.csdn.net/m0_56729804/article/details/130927733

阅读原文内容投诉

免责声明:

① 本站未注明“稿件来源”的信息均来自网络整理。其文字、图片和音视频稿件的所属权归原作者所有。本站收集整理出于非商业性的教育和科研之目的,并不意味着本站赞同其观点或证实其内容的真实性。仅作为临时的测试数据,供内部测试之用。本站并未授权任何人以任何方式主动获取本站任何信息。

② 本站未注明“稿件来源”的临时测试数据将在测试完成后最终做删除处理。有问题或投稿请发送至: 邮箱/279061341@qq.com QQ/279061341

软考中级精品资料免费领

  • 历年真题答案解析
  • 备考技巧名师总结
  • 高频考点精准押题
  • 2024年上半年信息系统项目管理师第二批次真题及答案解析(完整版)

    难度     813人已做
    查看
  • 【考后总结】2024年5月26日信息系统项目管理师第2批次考情分析

    难度     354人已做
    查看
  • 【考后总结】2024年5月25日信息系统项目管理师第1批次考情分析

    难度     318人已做
    查看
  • 2024年上半年软考高项第一、二批次真题考点汇总(完整版)

    难度     435人已做
    查看
  • 2024年上半年系统架构设计师考试综合知识真题

    难度     224人已做
    查看

相关文章

发现更多好内容

猜你喜欢

AI推送时光机
位置:首页-资讯-后端开发
咦!没有更多了?去看看其它编程学习网 内容吧
首页课程
资料下载
问答资讯