PythonCode：用于农机轨迹处理

import pandas as pda import numpy as npy import math as ma import time import datetime import os # 这两个是用于的打开文件选择框 from tkinter import Tk from tkinter.filedialog import askopenfilename # 这些是用来读写excel文件的 import xlwt import xlrd from xlutils.copy import copy ''' 2.0 版本相对 1.0 版本 更新的内容 : 1、用墨卡托坐标代替原来计算的XY坐标进行各种量的其他计算，plot画图所用的坐标还是 文献中简化后的 xy 坐标 2、标记了信号漂移点，并进行单独存储，结果文件中进行了删除处理 3、将文件进行了模块化处理 3.0 版本相对 2.0 版本 更新的 内容有： 1、将plot画图用的坐标换成 墨卡托坐标 ''' class GetPoint: def __init__(self, longtitude, latitude): self.latitude = latitude # 纬度 self.longtitude = longtitude # 经度 def ToY(self, weidu): a = 6378137 # 赤道圆的平均半径为 6378137 m b = 6356752 # 在WGS-84 中半短轴的长度为 6356752 m a1 = ma.pow(a, 2) b1 = ma.pow(b, 2) BO1 = ma.pow(ma.tan(self.latitude), 2) BO2 = ma.pow(1 / ma.tan(self.latitude), 2) Xo = a1 / ma.sqrt(b1 * BO1 + a1) # 参考点的横坐标 Yo = b1 / ma.sqrt(a1 * BO2 + b1) # 参考点的纵坐标 AO1 = ma.pow(ma.tan(weidu), 2) AO2 = ma.pow(1 / ma.tan(weidu), 2) Xa = a1 / ma.sqrt(b1 * AO1 + a1) # 目标点的横坐标 Ya = b1 / ma.sqrt(a1 * AO2 + b1) # 目标点的纵坐标 Y = ma.sqrt(ma.pow(Xa - Xo, 2) + ma.pow(Ya - Yo, 2)) # 目标点与参考点在平面坐标系中纵轴方向的距离 if weidu < self.latitude: Y = -Y return Y def ToX(self, jingdu): a = 6378137 b = 6356752 a1 = ma.pow(a, 2) b1 = ma.pow(b, 2) BO1 = ma.pow(ma.tan(self.latitude), 2) Xo = a1 / ma.sqrt(b1 * BO1 + a1) X = Xo * ma.fabs(jingdu - self.longtitude) # 目标点与参考点之间在平面坐标系中横轴方向上的距离 if jingdu < self.longtitude: X = -X return X def CalDisByLoLangitude(lon1, lat1, lon2, lat2): """ 根据两点的经纬度计算两点直线距离 :param lon1: 第一个点的 经度 :param lat1: 第一个点的 纬度 :param lon2: 第二个点的 经度 :param lat2: 第二个点的 纬度 :return: 两点的直线距离 """ distence = 0 temp = ma.sin(lat1) * ma.sin(lat2) * ma.cos(lon1 - lon2) + ma.cos(lat1) * ma.cos(lat2) if temp > 1 or temp < -1: # print('i = {},temp = {}'.format(i,temp)) if lon1 == lon2 and lat1 == lat2: distence = 0 else: distence = 6371004 * ma.acos(temp) * ma.pi / 180.0 return distence def LoLatitudeChangeToXY(origin, lon, lat): """ 一个经纬度坐标 转换 成 指定参考系中的 XY 坐标 :param origin: 原点，即参考经纬度 :param lon: 需要装换的 经度 :param lat: 需要装换的 纬度 :return: x y """ return origin.ToX(lon), origin.ToY(lat) def LoLatitudeChangeToXY_batch(b_lon, b_lat): """ 批量 经纬度坐标 转换 成 XY 坐标 :param b_lon: 经度 :param b_lat: 纬度 :return: x , y """ col = len(b_lon) # 数据长度 origin = GetPoint(b_lon[0], b_lat[0]) # 置 第一个 经纬度 为参考点 x = [] y = [] for i in range(1, col): x.append(origin.ToX(b_lon[i])) y.append(origin.ToY(b_lat[i])) return x, y def LonLatitude2WebMercator(lon, lat): """ 将经纬度转换成墨卡托坐标 :param lon: 经度 :param lat: 纬度 :return: 墨卡托坐标 X Y """ earthRad = 6378137.0 x = lon * npy.pi / 180 * earthRad a = lat * npy.pi / 180 y = earthRad / 2 * npy.log((1.0 + npy.sin(a)) / (1.0 - npy.sin(a))) return x, y def DisOnlineFromOnePointToAnother(lon1, lat1, angle1, lon2, lat2): """ 点2 到 点1 所在直线（直线斜率由 点1 平均航向角所得）的 最短距离 :param lon1: 点1 经度 :param lat1: 点1 纬度 :param angle1: 点1 平均航向角 :param lon2: 点2 经度 :param lat2: 点2 纬度 :return: 距离 """ origin = GetPoint(lon1, lat1) # 将 点1 设置为原点 x, y = LoLatitudeChangeToXY(origin, lon2, lat2) # 换算 点2 的 XY 大地坐标 # 求点2 的距离坐标 dis = CalDisByLoLangitude(lon1, lat1, lon2, lat2) # 计算 点1 到 点2 的直线距离 单位：米 # 点2 到 点1 连线的斜率 if x == 0: x_m = 0 y_m = dis else: k = y / x ang_k = ma.atan(k) x_m = dis * ma.cos(ang_k) # 单位 米 y_m = dis * ma.sin(ang_k) # 求点1 所在线 的斜率 if angle1 == 0 or angle1 == 360 and angle1 == 180: return ma.fabs(x) elif 0 < angle1 < 180: ag = (90 - angle1) * ma.pi / 180 # 倾斜角 弧度 else: ag = (270 - angle1) * ma.pi / 180 a = ma.tan(ag) # 斜率 # 直线方程 a X - Y = 0 return ma.fabs(a * x_m - y_m) / ma.sqrt(a * a + 1) def DisOnlineFromOnePointToAnotherByXY(lon1, lat1, angle1, lon2, lat2): """ 点2 到 点1 所在直线（直线斜率由 点1 平均航向角所得）的 最短距离 :param lon1: 点1 经度 :param lat1: 点1 纬度 :param angle1: 点1 平均航向角 :param lon2: 点2 经度 :param lat2: 点2 纬度 :return: 点到线的 最短距离 """ x1, y1 = LonLatitude2WebMercator(lon1, lat1) x2, y2 = LonLatitude2WebMercator(lon2, lat2) x = x2 - x1 y = y2 - y1 # 求点2 的距离坐标 dis = CalDisByLoLangitude(lon1, lat1, lon2, lat2) # 计算 点1 到 点2 的直线距离 单位：米 # 点2 到 点1 连线的斜率 if x == 0: x_m = 0 if y > 0: y_m = dis else: y_m = -dis else: k = y / x ang_k = ma.atan(k) x_m = dis * ma.cos(ang_k) # 单位 米 y_m = dis * ma.sin(ang_k) # 求点1 所在线 的斜率 if angle1 == 0 or angle1 == 360 and angle1 == 180: return ma.fabs(y_m) elif angle1 > 0 and angle1 < 180: ag = (90 - angle1) * ma.pi / 180 # 倾斜角 弧度 else: ag = (270 - angle1) * ma.pi / 180 a = ma.tan(ag) # 斜率 # 直线方程 a X - Y = 0 return ma.fabs(a * x_m - y_m) / ma.sqrt(a * a + 1) def GetSectionFileName(sourFileName, GpsTime, tpyeFlag=1): """ 获取切片文件名 根据原文件名，和 切片文件第一个GPS时间，以及类型标志（1=地块、0=道路，默认是地块）拼写切片文件名 :param sourFileName: 原文件名称 :param GpsTime: 切片文件第一个轨迹点的GPS时间 :param tpyeFlag: 类型标志:1-地块，0-道路 :return: 切片文件名 """ monthAndDay = "{:02d}{:02d}".format(GpsTime.month, GpsTime.day) hourAndMinute = "{:02d}{:02d}".format(GpsTime.hour, GpsTime.minute) sf = sourFileName.split('.') if tpyeFlag == 1: sectionFileName = sf[0] + '==' + monthAndDay + '-' + hourAndMinute + '-field.' + sf[1] else: sectionFileName = sf[0] + '==' + monthAndDay + '-' + hourAndMinute + '-road.' + sf[1] return sectionFileName def SplitfieldByMoreAtrributes(data, fileName): """ 联合多种属性分割地块 :param data: 读取的excel文件数据 :param fileName : excel 文件所在的路径及文件名 :return: data 处理后的excel文件数据 x 用于画图的横坐标 y 用于画图的纵坐标 """ head_text = list(data.columns) # 获取数据的 列标签 # 取列标签下标 # 标签名用字典表示 dict = {'depth': 0, 'toolState': 0, 'longitude': 0, 'latitude': 0, 'angel': 0} for i, tag in enumerate(head_text): if tag.find('作业深度') != -1: dict