EarthSatellite下的Skyfield可见区域

Question

我如何计算EarthSatellite下面的面积，以便在卫星经过时绘制覆盖的土地条带？

在Skyfield中有什么可以促进这一点的吗？

编辑:我只是想澄清一下我所说的卫星下面的区域是什么意思。考虑到地球是一个椭球体，我需要画出卫星下方可能观察到的最大面积。我知道如何绘制卫星路径，但现在我需要绘制一些线条来表示卫星在地球上空飞行时的可见区域。

Answer 1

你的编辑清楚地表明了你想要什么。卫星的可见面积可以很容易地计算出来(当地球被看作球体时)。可以在here上找到一些关于可见部分的背景信息。当地球被看作是扁球体时，计算可见面积将会困难得多(甚至可能是不可能的)。我认为更好的做法是修改这部分问题，并将其发布在数学上。

如果你想计算地球被看作球体时的可见面积，我们需要在Skyfield中做一些调整。通过使用TLE api加载卫星，您可以很容易地获得地球上的位置的子点。图书馆将其称为Geocentric位置，但实际上它是Geodetic位置(地球被视为扁球体)。为了纠正这个问题，我们需要调整Geocentric类的subpoint，使其使用Geocentric位置而不是Geodetic位置的计算。由于reverse_terra函数中的错误和缺少信息，我们还需要替换该函数。我们需要能够恢复地球的半径。这将导致以下结果：

from skyfield import api
from skyfield.positionlib import ICRF, Geocentric
from skyfield.constants import (AU_M, ERAD, DEG2RAD,
                                IERS_2010_INVERSE_EARTH_FLATTENING, tau)
from skyfield.units import Angle

from numpy import einsum, sqrt, arctan2, pi, cos, sin

def reverse_terra(xyz_au, gast, iterations=3):
    """Convert a geocentric (x,y,z) at time `t` to latitude and longitude.
    Returns a tuple of latitude, longitude, and elevation whose units
    are radians and meters.  Based on Dr. T.S. Kelso's quite helpful
    article "Orbital Coordinate Systems, Part III":
    https://www.celestrak.com/columns/v02n03/
    """
    x, y, z = xyz_au
    R = sqrt(x*x + y*y)

    lon = (arctan2(y, x) - 15 * DEG2RAD * gast - pi) % tau - pi
    lat = arctan2(z, R)

    a = ERAD / AU_M
    f = 1.0 / IERS_2010_INVERSE_EARTH_FLATTENING
    e2 = 2.0*f - f*f
    i = 0
    C = 1.0
    while i < iterations:
        i += 1
        C = 1.0 / sqrt(1.0 - e2 * (sin(lat) ** 2.0))
        lat = arctan2(z + a * C * e2 * sin(lat), R)
    elevation_m = ((R / cos(lat)) - a * C) * AU_M
    earth_R = (a*C)*AU_M
    return lat, lon, elevation_m, earth_R

def subpoint(self, iterations):
    """Return the latitude an longitude directly beneath this position.

    Returns a :class:`~skyfield.toposlib.Topos` whose ``longitude``
    and ``latitude`` are those of the point on the Earth's surface
    directly beneath this position (according to the center of the
    earth), and whose ``elevation`` is the height of this position
    above the Earth's center.
    """
    if self.center != 399:  # TODO: should an __init__() check this?
        raise ValueError("you can only ask for the geographic subpoint"
                            " of a position measured from Earth's center")
    t = self.t
    xyz_au = einsum('ij...,j...->i...', t.M, self.position.au)
    lat, lon, elevation_m, self.earth_R = reverse_terra(xyz_au, t.gast, iterations)

    from skyfield.toposlib import Topos
    return Topos(latitude=Angle(radians=lat),
                    longitude=Angle(radians=lon),
                    elevation_m=elevation_m)

def earth_radius(self):
    return self.earth_R

def satellite_visiable_area(earth_radius, satellite_elevation):
    """Returns the visible area from a satellite in square meters.

    Formula is in the form is 2piR^2h/R+h where:
        R = earth radius
        h = satellite elevation from center of earth
    """
    return ((2 * pi * ( earth_radius ** 2 ) * 
            ( earth_radius + satellite_elevation)) /
            (earth_radius + earth_radius + satellite_elevation))


stations_url = 'http://celestrak.com/NORAD/elements/stations.txt'
satellites = api.load.tle(stations_url)
satellite = satellites['ISS (ZARYA)']
print(satellite)

ts = api.load.timescale()
t = ts.now()

geocentric = satellite.at(t)
geocentric.subpoint = subpoint.__get__(geocentric, Geocentric)
geocentric.earth_radius = earth_radius.__get__(geocentric, Geocentric)

geodetic_sub = geocentric.subpoint(3)

print('Geodetic latitude:', geodetic_sub.latitude)
print('Geodetic longitude:', geodetic_sub.longitude)
print('Geodetic elevation (m)', int(geodetic_sub.elevation.m))
print('Geodetic earth radius (m)', int(geocentric.earth_radius()))

geocentric_sub = geocentric.subpoint(0)
print('Geocentric latitude:', geocentric_sub.latitude)
print('Geocentric longitude:', geocentric_sub.longitude)
print('Geocentric elevation (m)', int(geocentric_sub.elevation.m))
print('Geocentric earth radius (m)', int(geocentric.earth_radius()))
print('Visible area (m^2)', satellite_visiable_area(geocentric.earth_radius(), 
                                                    geocentric_sub.elevation.m))