根据经纬度坐标计算两点之间线段的交点

根据经纬度坐标计算两点之间线段的交点

class GisCheckUtils:
    L = 6381372 * math.pi * 2
    # 平面展开后,x轴等于周长
    W = L
    H = L / 2
    mill = 2.3

    def millierConvertion(self, lat, lon):
        # print(lat, lon)
        # 将经度从度数转换为弧度
        x = lon * math.pi / 180
        #  将纬度从度数转换为弧度
        y = lat * math.pi / 180
        # 米勒投影的转换
        # print(y)
        # print(math.log(math.tan(0.25 * math.pi + 0.4 * y)))
        y = 1.25 * math.log(math.tan(0.25 * math.pi + 0.4 * y))
        x = (self.W / 2) + (self.W / (2 * math.pi)) * x
        y = (self.H / 2) - (self.H / (2 * self.mill)) * y
        return [x, y]

    def xyToLatLon(self, x, y):
        x = (x - (self.W / 2)) / (self.W / (2 * math.pi))
        y = -1 * (y - (self.H / 2)) / (self.H / (2 * self.mill))

        y = (math.atan(math.pow(2.7182818284590452354, y / 1.25)) - 0.25 * math.pi) / 0.4
        lon = 180 / math.pi * x
        lat = 180 / math.pi * y
        return [lon, lat]

    # 排断两条直线是否相交
    def getIntersectPoint2(self, Points):
        p1 = Points[0]
        p2 = Points[1]
        p3 = Points[2]
        p4 = Points[3]
        A1 = p1[1] - p2[1]
        B1 = p2[0] - p1[0]
        C1 = A1 * p1[0] + B1 * p1[1]
        A2 = p3[1] - p4[1]
        B2 = p4[0] - p3[0]
        C2 = A2 * p3[0] + B2 * p3[1]
        det_k = A1 * B2 - A2 * B1

        if (math.fabs(det_k) < 0.00001):
            return None
        a = B2 / det_k
        b = -1 * B1 / det_k
        c = -1 * A2 / det_k
        d = A1 / det_k
        x = a * C1 + b * C2
        y = c * C1 + d * C2
        return [x, y]

    def segmentLatLonToPoint(self, line1, line2):
        points = []
        oneStartLat, oneStartLon, oneEndLat, oneEndLon = line1
        twoStartLat, twoStartLon, twoEndLat, twoEndLon = line2
        # print(twoStartLat, twoStartLon, twoEndLat, twoEndLon)
        points1 = self.millierConvertion(oneStartLon, oneStartLat)
        points2 = self.millierConvertion(oneEndLon, oneEndLat)
        points3 = self.millierConvertion(twoStartLon, twoStartLat)
        points4 = self.millierConvertion(twoEndLon, twoEndLat)
        points.append(points1)
        points.append(points2)
        points.append(points3)
        points.append(points4)
        return points

    def segIntersect(self, line1, line2):
        # 验证两条线有没有相交
        list = self.segmentLatLonToPoint(line1, line2)
        if self.segIntersect2(list) > 0:
            return True
        return False

    def getIntersectPoint(self, line1, line2):
        list = self.segIntersect(line1, line2)
        if list != True:
            return None
        # 获取两条直线相交的点
        points = self.segmentLatLonToPoint(line1, line2)
        latlon = self.getIntersectPoint2(points)
        if latlon:
            # return self.xyToLatLon(latlon[0], latlon[1])
            print(self.xyToLatLon(latlon[0], latlon[1]))
            return [latlon[0], latlon[1]]
        else:
            return None

    def segIntersect2(self, points):
        A = points[0]
        B = points[1]
        C = points[2]
        D = points[3]

        if (math.fabs(B[1] - A[1]) + math.fabs(B[0] - A[0]) + math.fabs(D[1] - C[1])
                + math.fabs(D[0] - C[0]) == 0):
            if (C[0] - A[0]) + (C[1] - A[1]) == 0:
                print("ABCD是同一个点!")
            else:
                print("AB是一个点,CD是一个点,且AC不同!")
            return 0

        if math.fabs(B[1] - A[1]) + math.fabs(B[0] - A[0]) == 0:
            if (A[0] - D[0]) * (C[1] - D[1]) - (A[1] - D[1]) * (C[0] - D[0]) == 0:
                print("A、B是一个点,且在CD线段上!")
            else:
                print("A、B是一个点,且不在CD线段上!")
            return 0

        if math.fabs(D[1] - C[1]) + math.fabs(D[0] - C[0]) == 0:
            if (D[0] - B[0]) * (A[1] - B[1]) - (D[1] - B[1]) * (A[0] - B[0]) == 0:
                print("C、D是一个点,且在AB线段上!")
            else:
                print("C、D是一个点,且不在AB线段上!")
            return 0

        if (B[1] - A[1]) * (C[0] - D[0]) - (B[0] - A[0]) * (C[1] - D[1]) == 0:
            print("线段平行,无交点x")
            return 0

        x = ((B[0] - A[0]) * (C[0] - D[0]) * (C[1] - A[1]) - C[0] * (B[0] - A[0]) * (C[1] - D[1]) + A[0] * (
                B[1] - A[1]) * (C[0] - D[0])) / ((B[1] - A[1]) * (C[0] - D[0]) - (B[0] - A[0]) * (C[1] - D[1]))
        y = ((B[1] - A[1]) * (C[1] - D[1]) * (C[0] - A[0]) - C[1] * (B[1] - A[1]) * (C[0] - D[0]) + A[1] * (
                B[0] - A[0]) * (C[1] - D[1])) / ((B[0] - A[0]) * (C[1] - D[1]) - (B[1] - A[1]) * (C[0] - D[0]))
        intersection = [x, y]
        if (intersection[0] - A[0]) * (intersection[0] - B[0]) <= 0 and (intersection[0] - C[0]) * (
                intersection[0] - D[0]) <= 0 and (intersection[1] - A[1]) * (intersection[1] - B[1]) <= 0 and (
                intersection[1] - C[1]) * (intersection[1] - D[1]) <= 0:
            if (A[0] == C[0] and A[1] == C[1]) or (A[0] == D[0] and A[1] == D[1]) or (
                    B[0] == C[0] and B[1] == C[1]) or (B[0] == D[0] and B[1] == D[1]):
                print("线段相交于端点上")
                return 2
            else:
                print("线段相交于点X:", intersection[0], " Y:", intersection[1])
                return 1
        else:
            # print("线段相交于虚交点(" + str(intersection[0]) + ","
            #       + str(intersection[1]) + ")!")
            return -1
            
gcu = GisCheckUtils()
line1 = oneStartLat, oneStartLon, oneEndLat, oneEndLon
line2 = twoStartLat, twoStartLon, twoEndLat, twoEndLon
# 交点坐标延长线不算
print(gcu.getIntersectPoint(line1, line2))
  • 0
    点赞
  • 3
    收藏
    觉得还不错? 一键收藏
  • 1
    评论
评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值