3

I'm trying to come up with a function where...
Input: Geodesic distance in miles or km
Output: The euclidean distance between any two gps points that are the input distance apart

I feel like I have some of the components

import numpy as np
from numpy import linalg as LA
from geopy.distance import geodesic

loc1 = np.array([40.099993, -83.166000])
loc2 = np.array([40.148652, -82.903962])

This is the euclidean distance between those two points

LA.norm(loc1-loc2)
#0.2665175636332336

This is the geodesic distance in miles between those two points

geodesic(loc1,loc2).miles
#14.27909749425243

My brain is running low on juice right now, anyone have any ideas on how I can make a function like:

geodesic_to_euclidean(14.27909749425243)
#0.2665175636332336
Jamalan
  • 482
  • 4
  • 15
  • 1
    Do you want the distance of the arc connecting two geolocations , assuming the earth's surface is a sphere? – Blake Oct 29 '20 at 03:35
  • I'm only looking at ranges within 50ish miles, so I think it would be fair to assume even flat surface? – Jamalan Oct 29 '20 at 03:37
  • `LA.norm` will not give you the correct value anyway because you're not working with cartesian coords. you need to calculate the angle between the two locations with respect to the center of the earth (in radians), then multiply by the radius of the earth. – Aaron Oct 29 '20 at 03:42
  • Well so I'm using the euclidean distance as a parameter within a clustering algorithm, DBSCAN. The only dimensions input into the algorithm are latitude and longitude. So I don't need the euclidean distance to represent anything real, just the 2-Norm distance between any two vectors. Is that problematic?...I'm super new to GIS stuff... – Jamalan Oct 29 '20 at 03:46
  • 1
    @Jamalan you'll loose accuracy the closer you get to the poles. This is the same reason the Mercator map projection blows up Antarctica much larger than it really is. – Aaron Oct 29 '20 at 03:48

1 Answers1

2

If you're okay with a great-circle distance, as mentioned in the comments, then this should work. It's the haversine distance:

def haversine(origin, destination, units='mi'):
    # Radian deltas
    origin_lat = radians(float(origin[0]))
    origin_lon = radians(float(origin[1]))
    destination_lat = radians(float(destination[0]))
    destination_lon = radians(float(destination[1]))
    lat_delta = destination_lat - origin_lat
    lon_delta = destination_lon - origin_lon

    # Radius of earth in meters
    r = 6378127

    # Haversine formula
    a = sin(lat_delta / 2) ** 2 + cos(origin_lat) * \
        cos(destination_lat) * sin(lon_delta / 2) ** 2
    c = 2 * asin(sqrt(a))
    meters_traveled = c * r

    scaling_factors = {
        "m:": 1,
        "km": 1 / 1000,
        "ft": 3.2808,  # meters to feet
        "mi:": 0.000621371  # meters to miles
    }

    return meters_traveled * scaling_factors[units]

If you already have the geodesic (great circle) distance in meters and you want the chord length, then you can do the following

def chord(geodesic_distance):
    """
    Chord length
    C = 2 * r * sin(theta/2)

    Arc length; which is geodesic distance in this case
    AL = R * theta

    therefore
    C = 2 * R * sin(AL/(2*R))
    """
    r = 6378127  # Radius of earth in meters

    return 2 * r * sin(geodesic_distance / (2 * r))
Blake
  • 986
  • 9
  • 24
  • 1
    for reference: [The haversine formula determines the great-circle distance between two points on a sphere given their longitudes and latitudes.](https://en.wikipedia.org/wiki/Haversine_formula) – Aaron Oct 29 '20 at 03:51
  • This is awesome thank you, definitely hanging on to this. Although this is haversine distance between two coordinates, vs geodesic distance as above. I was hoping the input could be just a distance like "14 miles". And it would return the euclidean distance between any two coordinates that were 14 miles apart. But based on Aaron's comment, I think I'm incorrectly assuming an evenly spaced grid. So there's not a global way to do this, excuse the pun – Jamalan Oct 29 '20 at 03:54
  • I see, so if you have the geodesic (arc curve) distance, you want to convert it to the straight-line distance? – Blake Oct 29 '20 at 03:59
  • 1
    @Jamalan definitely not an evenly spaced grid. I know it may not seem like what you're thinking about in your head, but what you're saying is to ask for the difference in length between a chord and an arc of the same angle. For small angles (small distances with respect to the radius of the earth) this difference will be minimal. – Aaron Oct 29 '20 at 03:59
  • Also, calculating the chord length from the lat-long coords isn't substantially easier than measuring the arc length using the haversine formula. This is again due to the fact that you don't have a square grid. – Aaron Oct 29 '20 at 04:02
  • @Blake, I want to convert it to the straight line distance in degrees... which sounds silly, but I'm trying to create a helper function that lets me experiment which making tighter and looser clusters. – Jamalan Oct 29 '20 at 04:03
  • 1
    @Jamalan `c` is that angle in radians already. just multiply by 180/pi – Aaron Oct 29 '20 at 04:04
  • @Aaron does any of this simplify if I'm only considering distances within a 50x50 mile area? – Jamalan Oct 29 '20 at 04:05
  • 1
    @Jamalan please see the above edit on chord length. Aaron is right though -- for nearby distances on Earth's surface, this will be minimally different – Blake Oct 29 '20 at 04:08
  • 1
    You guys are effing awesome. This is the good part of stack overflow. The non-cynical, just helpful people. Thank you both a ton. I will convert your wisdom and code into a sleek helper function. – Jamalan Oct 29 '20 at 04:10
  • @Jamalan I can't think of any optimizations other than throwing `numba` at the function. I would also speed test [`sklearn`](https://scikit-learn.org/stable/modules/generated/sklearn.metrics.pairwise.haversine_distances.html) and also try out some of the answers [here](https://stackoverflow.com/questions/29545704/fast-haversine-approximation-python-pandas) – Aaron Oct 29 '20 at 04:14