Мобильная версияИзначально мое распределение определяется N точками выборки и N соответствующими значениями скалярного поля. Поскольку точки расположены на сфере, в целом они расположены неравномерно. Это предполагает использование тепловой карты Вороного со следующей реализацией Scipy в сферическом случае.
Из объекта SphericalVoronoi Scipy можно получить доступ к списку регионов (ячейки Вороного) с декартовыми координатами каждой вершины ячейки. Их можно сопоставить с координатами длинной широты (используемыми проекцией Моллвейде в matplotlib), чтобы визуализировать регионы, используя код, близкий к приведенному ниже MWE.
Применение этого ко всем ячейкам приводит к путанице для регионов вблизи границы проекции Моллвейде. Эти проблемные регионы характеризуются наличием как минимум двух вершин с долготой противоположного знака (но близкой к pi по абсолютной величине, при условии, что N достаточно велико и точки выборки распределены равномерно). Вот пример, где такие регионы обозначаются красным крестом в точке генератора:
Чтобы решить эту проблему, моей первой идеей было бы сделать следующее:
- Идентифицировать все проблемные области (уже сделано, критерий с долготой)
- вдоль каждого края ячейки между двумя вершинами противоположной долготы добавьте две вершины с долготой np.pi + e и np.pi - e, где e = 1e-10 — некоторое небольшое число.
- Разделите каждую ячейку на два подмножества вершин с долготами противоположного знака, включая четыре, добавленные в предыдущем шаг
- Постройте каждое подмножество как одну ячейку Вороного, хотя оба подмножества соответствуют одному и тому же значению скалярного распределения, которое мы выбираем.
Есть ли у вас какие-либо предложения по упрощению процесса с учетом того, что я работаю с Python? Будем рады любым советам или идеям по улучшению!
Минимальный рабочий пример
Сначала мы импортируем необходимые пакеты и определяем нашу выборку в сфере:
import numpy as np
import matplotlib.pyplot as plt
from scipy.spatial import SphericalVoronoi
def long_lat_from_xyz(x, y, z):
lat = np.asin(z)
long = np.atan2(y, x)
return long, lat
def sample_points(nb_points):
# Sample points on the sphere and return their cartesian and long-lat coordinates
long_lat = np.zeros(shape=(nb_points, 2))
xyz = np.zeros(shape=(nb_points, 3))
ratio = (1 + np.sqrt(5))*np.pi/2
for i, x in enumerate(np.linspace(-1, 1, nb_points)):
r = np.sqrt(1 - x**2)
angle = i*ratio
y = r*np.cos(angle)
z = r*np.sin(angle)
x, y, z = y, z, x # Cosmetic change to have the z-axis as the vertical one
xyz[i, :] = [x, y, z]
long_lat[i, :] = long_lat_from_xyz(x, y, z)
return xyz, long_lat
Мы реализуем, как идентифицируются регионы на границе (в соответствии с тем, что предшествует) и как отображать регионы, здесь в простой схеме:
def are_points_on_opposite_sides(xyz1, xyz2):
x1, y1, z1 = xyz1
x2, y2, z2 = xyz2
long1, _ = long_lat_from_xyz(x1, y1, z1)
long2, _ = long_lat_from_xyz(x2, y2, z2)
return (long1*long2 < 0) and (np.abs(long1) > np.pi/2) and (np.abs(long2) > np.pi/2)
def is_region_at_boundary(region):
n = len(region)
for i in range(n):
v1 = sv.vertices[region]
v2 = sv.vertices[region][(i + 1) % n]
if are_points_on_opposite_sides(v1, v2):
return True
return False
def plot_region(region, ax):
for i in range(len(region)):
x_s, y_s, z_s = sv.vertices[region]
x_e, y_e, z_e = sv.vertices[region][(i + 1) % len(region)]
long_s, lat_s = long_lat_from_xyz(x_s, y_s, z_s)
long_e, lat_e = long_lat_from_xyz(x_e, y_e, z_e)
ax.plot([long_s, long_e], [lat_s, lat_e], color='k')
В конечном итоге мы построим всю диаграмму Вороного:
nb_points = 50
xyz, long_lat = sample_points(nb_points)
sv = SphericalVoronoi(xyz, 1, np.array([0, 0, 0])) # xyz, radius, center
sv.sort_vertices_of_regions()
fig = plt.figure()
ax = fig.add_subplot(111, projection='mollweide')
for n, region in enumerate(sv.regions):
if not is_region_at_boundary(region):
plot_region(region, ax)
else:
x, y, z = sv.points[n, :]
long, lat = long_lat_from_xyz(x, y, z)
ax.scatter(long, lat, marker='x', color='red')
plt.show()
Подробнее здесь: https://stackoverflow.com/questions/798 ... projection
