Как построить поле обзора спутника наблюдения Земли, когда он находится близко к полюсам, используя базовую карту?

Question

Я пытаюсь нарисовать максимальное (теоретическое) поле зрения спутника вдоль его орбиты. Я использую Базовую карту, на которой я хочу построить различные позиции вдоль орбиты (с разбросом), и я хотел бы добавить все поле зрения, используя метод тканей (или эквивалент). Приведенный ниже код работает нормально до тех пор, пока широта не достигнет примерно 75 градусов северной широты на высоте 300 км. За пределами которого код выводит ValueError: «ValueError: неопределенная обратная геодезическая (может быть антиподальной точкой)»

import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
import math

earth_radius = 6371000.  # m
fig = plt.figure(figsize=(8, 6), edgecolor='w')
m = Basemap(projection='cyl', resolution='l',
            llcrnrlat=-90, urcrnrlat=90,
            llcrnrlon=-180, urcrnrlon=180)

# draw the coastlines on the empty map
m.drawcoastlines(color='k')

# define the position of the satellite
position = [300000., 75., 0.]  # altitude, latitude, longitude

# radius needed by the tissot method
radius = math.degrees(math.acos(earth_radius / (earth_radius + position[0])))
m.tissot(position[2], position[1], radius, 100, facecolor='tab:blue', alpha=0.3)
m.scatter(position[2], position[1], marker='*', c='tab:red')

plt.show()

Следует отметить, что код прекрасно работает на южном полюсе (широта ниже -75). Я знаю, что это известная ошибка, есть ли известное решение этой проблемы? Спасибо за вашу помощь!

Answer 1

То, что вы обнаружили, является некоторыми из ограничений Базовой карты. Давайте переключимся на Cartopy пока. Рабочий код будет другим, но не очень.

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import math

earth_radius = 6371000.
position = [300000., 75., 0.]   # altitude (m), lat, long
radius = math.degrees(math.acos(earth_radius / (earth_radius + position[0])))
print(radius)  # in subtended degrees??

fig = plt.figure(figsize=(12,8))

img_extent = [-180, 180, -90, 90]

# here, cartopy's' `PlateCarree` is equivalent with Basemap's `cyl` you use
ax = fig.add_subplot(1, 1, 1, projection = ccrs.PlateCarree(), extent = img_extent)

# for demo purposes, ...
# let's take 1 subtended degree = 112 km on earth surface (*** you set the value as needed ***)
ax.tissot(rad_km=radius*112, lons=position[2], lats=position[1], n_samples=64, \
             facecolor='red', edgecolor='black', linewidth=0.15, alpha = 0.3)

ax.coastlines(linewidth=0.15)
ax.gridlines(draw_labels=False, linewidth=1, color='blue', alpha=0.3, linestyle='--')
plt.show()

С кодом, приведенным выше, результирующий график:

Теперь, если мы используем ортографическую проекцию, (замените соответствующую строку кода этим)

ax = fig.add_subplot(1, 1, 1, projection = ccrs.Orthographic(central_longitude=0.0, central_latitude=60.0))

Вы получаете этот участок: