代码之家  ›  专栏  ›  技术社区  ›  jovicbg

如何在pandas数据帧python中找到gps坐标之间的夹角

  •  1
  • jovicbg  · 技术社区  · 6 年前

    我有带测量坐标和单元坐标的数据框。

    我需要找到连接这两点和北极的直线之间的每一行角度(方位角)。

    数据框:

    id     cell_lat     cell_long     meas_lat     meas_long
    1      53.543643    11.636235     53.44758     11.03720
    2      52.988823    10.0421645    53.03501     9.04165
    3      54.013442    9.100981      53.90384     10.62370
    

    我在网上找到了一些代码,但如果这真的能帮助我更接近解决方案的话,那就没有了。

    我用过 this 函数,但不确定是否正确,我想有更简单的解决方案。

    欢迎任何帮助或提示,提前谢谢。

    1 回复  |  直到 6 年前
        1
  •  1
  •   liamvharris    6 年前

    这个问题最棘手的部分是将大地坐标(纬度,经度)转换为笛卡尔坐标(x,y,z)。如果你看看 https://en.wikipedia.org/wiki/Geographic_coordinate_conversion 您可以看到如何做到这一点,其中包括选择一个参考系统。假设我们选择ECEF( https://en.wikipedia.org/wiki/ECEF ,以下代码将计算您要查找的角度:

    def vector_calc(lat, long, ht):
        '''
        Calculates the vector from a specified point on the Earth's surface to the North Pole.
        '''
        a = 6378137.0  # Equatorial radius of the Earth
        b = 6356752.314245  # Polar radius of the Earth
    
        e_squared = 1 - ((b ** 2) / (a ** 2))  # e is the eccentricity of the Earth
        n_phi = a / (np.sqrt(1 - (e_squared * (np.sin(lat) ** 2))))
    
        x = (n_phi + ht) * np.cos(lat) * np.cos(long)
        y = (n_phi + ht) * np.cos(lat) * np.sin(long)
        z = ((((b ** 2) / (a ** 2)) * n_phi) + ht) * np.sin(lat)
    
        x_npole = 0.0
        y_npole = 6378137.0
        z_npole = 0.0
    
        v = ((x_npole - x), (y_npole - y), (z_npole - z))
    
        return v
    
    def angle_calc(lat1, long1, lat2, long2, ht1=0, ht2=0):
        '''
        Calculates the angle between the vectors from 2 points to the North Pole.
        '''
        # Convert from degrees to radians
        lat1_rad = (lat1 / 180) * np.pi
        long1_rad = (long1 / 180) * np.pi
        lat2_rad = (lat2 / 180) * np.pi
        long2_rad = (long2 / 180) * np.pi
    
        v1 = vector_calc(lat1_rad, long1_rad, ht1)
        v2 = vector_calc(lat2_rad, long2_rad, ht2)
    
        # The angle between two vectors, vect1 and vect2 is given by:
        # arccos[vect1.vect2 / |vect1||vect2|]
        dot = np.dot(v1, v2)  # The dot product of the two vectors
        v1_mag = np.linalg.norm(v1)  # The magnitude of the vector v1
        v2_mag = np.linalg.norm(v2)  # The magnitude of the vector v2
    
        theta_rad = np.arccos(dot / (v1_mag * v2_mag))
        # Convert radians back to degrees
        theta = (theta_rad / np.pi) * 180
    
        return theta
    
    angles = []
    for row in range(df.shape[0]):
        cell_lat = df.iloc[row]['cell_lat']
        cell_long = df.iloc[row]['cell_long']
        meas_lat = df.iloc[row]['meas_lat']
        meas_long = df.iloc[row]['meas_long']
    
        angle = angle_calc(cell_lat, cell_long, meas_lat, meas_long)
    
        angles.append(angle)
    

    这将读取数据框中的每一行,计算角度并将其附加到列表角度。显然,在计算出这些角度后,你可以用它们做你喜欢的事情。

    希望能有帮助!