Python: расчет средней широты, включая пересечение линий дат

У меня есть массив (lons) значений долготы в диапазоне [-180, 180]. Мне нужно найти среднее значение временного ряда. Это легко сделать с

np.mean(lons)

Это прямое среднее значение, конечно, не работает, если ряд содержит значения по обе стороны от линии дат. Каков правильный способ вычисления среднего для всех возможных случаев? Обратите внимание: я бы предпочел не иметь условия, которое по-разному обрабатывает случаи пересечения линии дат.

Я поиграл с np.unwrap после преобразования градусов в рад, но я знаю, что мои расчеты неверны, потому что небольшой процент случаев дает мне средние долготы где-то около 0 градусов (меридиан) над Африкой. Это невозможно, так как это набор данных об океане.

Спасибо.

РЕДАКТИРОВАТЬ: Теперь я понимаю, что более точный способ вычисления среднего [широта, долгота] положения временного ряда может заключаться в преобразовании в декартову сетку. Я могу пойти по этому пути.


person InitialConditions    schedule 21.11.2017    source источник


Ответы (2)


Это приложение для направленной статистики, где среднее угловое значение вычисляется в комплексной плоскости (см. этого раздела). Результатом является комплексное число, мнимая часть которого представляет собой средний угол:

import numpy as np

def angular_mean(angles_deg):
    N = len(angles_deg)
    mean_c = 1.0 / N * np.sum(np.exp(1j * angles_deg * np.pi/180.0))
    return np.angle(mean_c, deg=True)

lons = [
    np.array([-175, -170, 170, 175]),  # broad distribution
    np.random.rand(1000)               # narrow distribution
]

for lon in lons:
    print angular_mean(lon), np.mean(lon)

Как видите, среднее арифметическое и среднее угловое очень похожи для узкого распределения, тогда как они значительно различаются для широкого распределения.

Использование декартовых координат нецелесообразно, так как центр масс будет находиться внутри Земли, но, поскольку вы используете данные о поверхности, я предполагаю, что вы хотите, чтобы он располагался на поверхности.

person sfinkens    schedule 21.11.2017
comment
Спасибо за ваше решение. Мне действительно удалось успешно решить это, используя преобразование в декартовы координаты, чтобы найти среднюю широту и долготу. Мне не нужно значение z, так что на самом деле это не трехмерный центр масс, а просто среднее значение x и y. Я проверил ваш код, и результаты идентичны или находятся в пределах 1% - person InitialConditions; 21.11.2017
comment
Я ошибся насчет трансформации. Долготы в [-180, 180] просто прекрасны - person sfinkens; 21.11.2017
comment
Но работал ли исходный код? Это было для меня. Вы просто говорите, что преобразование в [0, 360] и обратно было ненужным? - person InitialConditions; 21.11.2017
comment
Я думаю, что да - person sfinkens; 21.11.2017
comment
Что касается декартового метода: вы не можете восстановить широту без компонента z. широта ~ arcsin(z/R) - person sfinkens; 21.11.2017
comment
Нет, ты не можешь. Но я все еще могу получить то, что хочу, то есть среднюю широту и долготу набора координат [широта, долгота]. См. второй ответ ниже: stackoverflow.com/questions/1185408/ - person InitialConditions; 21.11.2017

Вот мое решение. Обратите внимание, что я вычисляю среднюю широту и долготу, а также среднее расстояние (mean_dist) координат [широта, долгота] из рассчитанной средней широты (lat_mean) и средней долготы (lon_mean). Причина в том, что меня также интересует, насколько сильно отличается центральная [широта, долгота]. Я считаю, что это правильно, но я открыт для обсуждения!

lat_size = np.size(lats)
lon_rad = np.deg2rad(lons)  # lons in degrees [-180, 180]
lat_rad = np.deg2rad(lats)  # lats in degrees [-90, 90]
R = 6371  #  Approx radius of Earth (km)    
x = R * np.cos(lat_rad) * np.cos(lon_rad)
y = R * np.cos(lat_rad) * np.sin(lon_rad)
z = R * np.sin(lat_rad)

x_mean = np.mean(x)
y_mean = np.mean(y)
z_mean = np.mean(z)    

lat_mean = np.rad2deg(np.arcsin(z_mean / R))
lon_mean = np.rad2deg(np.arctan2(y_mean, x_mean))

# Calculate distance from centre point for each [lat, lon] pair    
dist_list = np.empty(lat_size)
dist_list.fill(np.nan)
p = 0
for lat, lon in zip(lats, lons):
    coords_1 = (lat, lon)
    coords_2 = (lat_mean, lon_mean )
    dist_list[p] = geopy.distance.vincenty(coords_1, coords_2).km
    p = p + 1
mean_dist = np.mean(dist_list)
return lat_mean, lon_mean, mean_dist
person InitialConditions    schedule 21.11.2017