У меня есть данные в формате:
Az Elv Value
0 10 15
0 30 30
0 60 22
15 10 16
15 30 20
и т. Д.
Я пытаюсь интерполировать эти измерения, чтобы получить плавный контур неба и измерения. Мои данные не имеют измерений непосредственно над головой, а также измерений ниже горизонта.
В настоящее время я выполняю следующие действия:
- Преобразовать в декартову
- Создание диапазонов из минимальных и максимальных декартовых координат
- Сделать сетку используя
np.meshgrid
- Использование griddata для получения значений сетки
- Преобразовать координаты сетки обратно в полярные
- Построить контур, используя сетку тета, r и значения
Я получаю хороший полярный график интерполированных данных, которыми я доволен, за исключением заполненных значений в зените (от 90 градусов до 75 градусов). У меня там нет измерений, и я хочу, чтобы они оставались пустыми. Я не могу понять, как это сделать. Если вы посмотрите на график, то увидите, что интерполированные измерения в зените не имеют смысла.
Бонусный вопрос, я переворачиваю ось r так, что 90 находится в центре, а 0 - снаружи. Чтобы это выглядело как ночное небо. Могу ли я просто вычесть значения r из 90, чтобы получить это значение?
Ниже мой код:
data = ascii.read(file)
azimuth = np.array(data['az'])
elv = np.array(data['elv'])
values = np.array(data['nsb'])
theta = np.radians(azimuth)
r = elv
z = values
x = r * np.cos(theta)
y = r * np.sin(theta)
xgrid = np.linspace(x.min(), x.max(), len(x))
ygrid = np.linspace(y.min(), y.max(), len(y))
xgrid, ygrid = np.meshgrid(xgrid, ygrid)
zgrid = griddata((x, y), z, (xgrid, ygrid), method='cubic')
#zgrid = ma.masked_invalid(zgrid)
print(zgrid.min(), zgrid.max())
r = np.sqrt(xgrid**2 + ygrid**2)
#r = ma.masked_where(r < 15, r)
theta = np.arctan2(ygrid, xgrid)
ax = plt.subplot(111, projection='polar')
ax.set_facecolor('k')
colors = plt.cm.get_cmap('cubehelix_r')
levels, steps = np.linspace(18.375, 22.375, 35, retstep=True)
print(steps)
cax = ax.contourf(theta, 90-r, zgrid, levels=levels, cmap=colors)
cbar = plt.colorbar(cax)
cbar.set_label(r'mag/arcsec$^2$')
ax.set_theta_zero_location('N')
ax.set_theta_direction(-1)
ax.set_yticklabels('')
ax.yaxis.grid(False)
ax.xaxis.grid(False)
plt.ylim(0,90)
plt.show()
Спасибо!