14

Is there's a library or a way to calculate the center point for several geolocations points? This is my list of geolocations based in New York and want to find the approximate midpoint geolocation

L = [
     (-74.2813611,40.8752222),
     (-73.4134167,40.7287778),
     (-74.3145014,40.9475244),
     (-74.2445833,40.6174444),
     (-74.4148889,40.7993333),
     (-73.7789256,40.6397511)
    ]
4
  • 1
    Check out this question on the GIS stack exchange: gis.stackexchange.com/questions/12120/…
    – PFlans
    Jun 17, 2016 at 15:59
  • 1
    how about calculating the simple average of all coordinates (they are quite close to each other)?
    – gdlmx
    Jun 17, 2016 at 16:00
  • @PFlans that is a very helpful source! thank you!
    – mongotop
    Jun 17, 2016 at 16:07
  • @gdlmx yes, they are close. I will calculate the averages for lat and long as center point. thank you!
    – mongotop
    Jun 17, 2016 at 16:08

5 Answers 5

23

After the comments I received and comment from HERE

With coordinates that close to each other, you can treat the Earth as being locally flat and simply find the centroid as though they were planar coordinates. Then you would simply take the average of the latitudes and the average of the longitudes to find the latitude and longitude of the centroid.

lat = []
long = []
for l in L :
  lat.append(l[0])
  long.append(l[1])

sum(lat)/len(lat)
sum(long)/len(long)

-74.07461283333332, 40.76800886666667
0
9

Based on: https://gist.github.com/tlhunter/0ea604b77775b3e7d7d25ea0f70a23eb

Assume you have a pandas DataFrame with latitude and longitude columns, the next code will return a dictionary with the mean coordinates.

import math

x = 0.0
y = 0.0
z = 0.0

for i, coord in coords_df.iterrows():
    latitude = math.radians(coord.latitude)
    longitude = math.radians(coord.longitude)

    x += math.cos(latitude) * math.cos(longitude)
    y += math.cos(latitude) * math.sin(longitude)
    z += math.sin(latitude)

total = len(coords_df)

x = x / total
y = y / total
z = z / total

central_longitude = math.atan2(y, x)
central_square_root = math.sqrt(x * x + y * y)
central_latitude = math.atan2(z, central_square_root)

mean_location = {
    'latitude': math.degrees(central_latitude),
    'longitude': math.degrees(central_longitude)
    }
3

Considering that you are using signed degrees format (more), simple averaging of latitude and longitudes would create problems for even small regions near to antimeridian (i.e. + or - 180-degree longitude) due to discontinuity of longitude value at this line (sudden jump between -180 to 180).

Consider two locations whose longitudes are -179 and 179, their mean would be 0, which is wrong.

3

This link can be useful, first convert lat/lon into an n-vector, then find average. A first stab at converting the code into Python is below

import numpy as np
import numpy.linalg as lin

E = np.array([[0, 0, 1],
              [0, 1, 0],
              [-1, 0, 0]])

def lat_long2n_E(latitude,longitude):
    res = [np.sin(np.deg2rad(latitude)),
           np.sin(np.deg2rad(longitude)) * np.cos(np.deg2rad(latitude)),
           -np.cos(np.deg2rad(longitude)) * np.cos(np.deg2rad(latitude))]
    return np.dot(E.T,np.array(res))

def n_E2lat_long(n_E):
    n_E = np.dot(E, n_E)
    longitude=np.arctan2(n_E[1],-n_E[2]);
    equatorial_component = np.sqrt(n_E[1]**2 + n_E[2]**2 );
    latitude=np.arctan2(n_E[0],equatorial_component);
    return np.rad2deg(latitude), np.rad2deg(longitude)

def average(coords):
    res = []
    for lat,lon in coords:
        res.append(lat_long2n_E(lat,lon))
    res = np.array(res)
    m = np.mean(res,axis=0)
    m = m / lin.norm(m)
    return n_E2lat_long(m)
        

n = lat_long2n_E(30,20)
print (n)
print (n_E2lat_long(np.array(n)))

# find middle of france and libya
coords = [[30,20],[47,3]]
m = average(coords)
print (m)

enter image description here

1
  • 1
    Thank you @BBSysDyn. This is a wonderful explanation!
    – mongotop
    Jul 7, 2021 at 1:52
1

I would like to improve on the @BBSysDyn'S answer. The average calculation can be biased if you are calculating the center of a polygon with extra vertices on one side. Therefore the average function can be replaced with centroid calculation explained here

def get_centroid(points):

    x = points[:,0]
    y = points[:,1]
    
    # Solving for polygon signed area
    A = 0
    for i, value in enumerate(x):
        if i + 1 == len(x):
            A += (x[i]*y[0] - x[0]*y[i])
        else:
            A += (x[i]*y[i+1] - x[i+1]*y[i])
    A = A/2

    #solving x of centroid
    Cx = 0
    for i, value in enumerate(x):
        if i + 1 == len(x):
            Cx += (x[i]+x[0]) * ( (x[i]*y[0]) - (x[0]*y[i]) )
        else:
            Cx += (x[i]+x[i+1]) * ( (x[i]*y[i+1]) - (x[i+1]*y[i]) )
    Cx = Cx/(6*A)

    #solving y of centroid
    Cy = 0
    for i , value in enumerate(y):
        if i+1 == len(x):
            Cy += (y[i]+y[0]) * ( (x[i]*y[0]) - (x[0]*y[i]) )
        else:
            Cy += (y[i]+y[i+1]) * ( (x[i]*y[i+1]) - (x[i+1]*y[i]) )
    Cy = Cy/(6*A)

    return Cx, Cy

Note: If it is a polygon or more than 2 points, they must be listed in order that the polygon or shape would be drawn.

1
  • Thanks! it works really well :), for vectors with 2 or 3 values A sometimes is equal to 0 so it breaks, I added a flag in my implementation to check that.
    – ruloweb
    Jul 1 at 0:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.