Search code examples
pythondistancebearingazimuth

Calculate Harvesine and initial bearing for successive points of a panda DataFrame in python


There are so many packages that provide this calculation although most of them are based on points and not data frames or maybe I am making a mistake! I have found this method that workers with my panda dataframe of Latitude and Longitude columns:

def haversine(lat1, lon1, lat2, lon2, to_radians=True, earth_radius=6378137):
   """
   slightly modified version: of http://stackoverflow.com/a/29546836/2901002

   Calculate the great circle distance between two points
   on the earth (specified in decimal degrees or in radians)

   All (lat, lon) coordinates must have numeric dtypes and be of equal length.
   """
   if to_radians:
       lat1, lon1, lat2, lon2 = map(np.radians, [lat1, lon1, lat2, lon2])
       a = np.sin((lat2-lat1)/2.0)**2 + \
           np.cos(lat1) * np.cos(lat2) * np.sin((lon2-lon1)/2.0)**2
   return earth_radius * 2 * np.arcsin(np.sqrt(a))

but all initial bearing or azimuths I tried, don't accept the dataframe series and trying the numpy array will still return me zeros! Is there a certain way to do so for a successive rows of a dataframe? I would like to calculate initial bearing between successive points. In R the bearing function will do the job with a dataframe just wondering if there is an equivalent in Python.


Solution

  • Update: I found the problem. I was using the R method to be able to find bearing between successive rows so I was basically removing the first row and the last row making two sets of dataframes with two columns but it worked perfectly with shift() and I wrote my own bearing function which was easier than using the one out there... So I made the tow dataframes below from pts my main dataframe: latlon_a = pts latlon_b = pts.shift() and my own initial bearing function:

    def initial_bearing(lon1, lat1, lon2, lat2):
       """
       My own version based on R source
    
       Calculate the initial bearing between two points
    
       All (latitude, longitude) coordinates must have numeric dtypes and be of equal length.
       """
       lat1, lon1, lat2, lon2 = map(np.radians, [lon1, lat1, lon2, lat2])
       delta1 = lon1-lon2
       term1 = np.sin(delta1) * np.cos(lat2)
       term2 = np.cos(lat1) * np.sin(lat2)
       term3 = np.sin(lat1) * np.cos(lat2) * np.cos(delta1)
       rad = np.arctan2(term1, (term2-term3))
       bearing = np.rad2deg(rad)
       return (bearing + 360) % 360
    
    
    bearing = initial_bearing(latlon_a['longitude'],latlon_a['latitude'],
                              latlon_b['longitude'],latlon_b['latitude'])
    

    This worked perfectly for me and results back the initial bearing. for funial bearing you can just replace or add the line below to return: return (bearing + 180) % 360