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

Я пытаюсь нарисовать максимальное (теоретическое) поле зрения спутника по его орбите. Я использую Basemap, на котором я хочу нанести разные положения по орбите (с разбросом), и я хотел бы добавить все поле зрения, используя метод tissot (или аналог). Приведенный ниже код работает нормально, пока широта не достигнет примерно 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). Я знаю, что это известная ошибка, есть ли известное решение этой проблемы? Спасибо за вашу помощь!


python matplotlib-basemap cartopy satellite fieldofview
person Flo    schedule 15.01.2019    source источник

Ответы (1)


arrow_upward
3
arrow_downward

Вы обнаружили некоторые ограничения Basemap. Давайте пока переключимся на 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))

вы получаете этот сюжет:

введите здесь описание изображения

person swatchai    schedule 16.01.2019
comment
Спасибо за ваш ответ! Я считаю, что это то, чего я хочу. Попробую Картопи. Ваше здоровье - person Flo; 16.01.2019
comment
Я немного опоздал к этому обсуждению, но переменный радиус — это угол, поэтому alpha = math.acos(earth_radius / (earth_radius + position[0])); радиус = альфа * земной_радиус; так что сюжет не должен быть таким: ax.tissot(rad_km=radius/1000.....? - person rul30; 22.10.2020