import logging logging.basicConfig(level=logging.INFO, format='%(asctime)s-%(filename)s[line:%(lineno)d]-%(levelname)s:%(message)s', datefmt='%Y-%m-%d %H:%M:%S') import numpy as np from shapely.geometry import LineString, Point from shapely.ops import substring def pnt2polar(p): # 平面点转极坐标 x, y = p pi = np.pi alpha = np.arctan(abs(y / x)) if x != 0. else 0. if x > 0: if y > 0: angle = alpha elif y == 0: angle = 0 else: angle = pi * 2 - alpha elif x == 0: if y > 0: angle = pi / 2 elif y == 0: angle = 0 else: angle = 3 * pi / 2 else: if y > 0: angle = pi - alpha elif y == 0: angle = pi else: angle = pi + alpha return round(angle, 3) def vline(line): # line shapely.geometry.LineString if line.length < 1e-6: res = None # todo 处理空线的情况 else: x0, y0 = line.coords[0] x1, y1 = line.coords[-1] vec = [x1-x0, y1-y0] res = pnt2polar(vec) return res def lineXroad(line=None, road=None, near=30, N=11): # 在不完全重合的情况下 判断一段线是否是另一条线的一部分 # line / road LineString # near 缓冲区半径 default 30m # N 参考点的数量 seg_pnts = [substring(road, _, _, normalized=1) for _ in np.linspace(0, 1, N)] # 取点 seg_line = [_.buffer(near).intersection(road) for _ in seg_pnts] # 截局部 seg_line = [_ if _.geom_type=='LineString' else _.geoms[0] for _ in seg_line] # 截出多线取第一部分 iline = [_.buffer(near).intersection(line) for _ in seg_pnts] # 截取线路 iline = [_ if _.geom_type=='LineString' else _.geoms[0] for _ in iline] rvec = [vline(_) for _ in seg_line] # road方向 逐段的 lvec = [vline(_) for _ in iline] # line方向 logging.warning(rvec) logging.warning(lvec) flag = [abs(m-n) if n is not None else -1 for m, n in zip(rvec, lvec)] res = sum([_ < 0.262 or _ > 6.02 if _ >= 0 else 0 for _ in flag]) return res > N//2+1 if __name__ == "__main__": line = [(0, 0), (3, 0), (3, 3), (6, 3), (6, 0), (9, 0)] road1 = [(-1, -0.1), (2, -0.1)] # True 超过一半重合 road2 = [(2.95, 0), (2.95, 3)] # True 完全重合 road3 = [(6, 3.05), (3, 3.05)] # False 方向相反 road4 = [(6, -0.1), (9, -0.1), (15, 0)] # false 少部分重合 roads = [LineString(_) for _ in [road1, road2, road3, road4]] line = LineString(line) for road in roads: res = lineXroad(line=line, road=road, near=0.3) logging.warning(res)
