Skip to content

Instantly share code, notes, and snippets.

@kerel-fs

kerel-fs/iss_transit_astropy.py

Forked from cbassa/iss_transit_astropy.py
Last active November 26, 2019 03:05
Show Gist options
  • Save kerel-fs/0a0f40cc7e693dadd075b0ae48fdcd28 to your computer and use it in GitHub Desktop.
Save kerel-fs/0a0f40cc7e693dadd075b0ae48fdcd28 to your computer and use it in GitHub Desktop.
Download ZIP
Reproducing the ISS->Venus->Solar transit of June 8, 2004 (Astropy & skyfueld & Astropy+Skyfield Frame Trafos)
Raw
iss_transit_astropy.py
#!/usr/bin/env python3
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Circle
import astropy.units as u
from astropy.time import Time
from astropy.coordinates import get_body, solar_system_ephemeris
from astropy.coordinates import EarthLocation, GCRS, FK5
from astropy.coordinates import SkyCoord, PrecessedGeocentric, CartesianRepresentation
from astropy.wcs import wcs
from astropy.coordinates.earth_orientation import precession_matrix_Capitaine, nutation_components2000B
from astropy.coordinates.matrix_utilities import rotation_matrix, matrix_product, matrix_transpose
from astropy.coordinates import cartesian_to_spherical
from sgp4.earth_gravity import wgs84
from sgp4.io import twoline2rv
# From astropy, but with units
def nutation_matrix(epoch):
"""
Nutation matrix generated from nutation components.
Matrix converts from mean coordinate to true coordinate as
r_true = M * r_mean
"""
# TODO: implement higher precision 2006/2000A model if requested/needed
epsa, dpsi, deps = nutation_components2000B(epoch.jd) # all in radians
epsa *= u.rad
dpsi *= u.rad
deps *= u.rad
return matrix_product(rotation_matrix(-(epsa + deps), 'x'),
rotation_matrix(-dpsi, 'z'),
rotation_matrix(epsa, 'x'))
if __name__ == "__main__":
# Set time
t = Time("2004-06-08T10:09:17", scale="utc", format="isot")
dts = np.arange(5)*u.s
# Generate timestamps
tms = list([t+dt for dt in dts])
# Set location
loc = EarthLocation(lat=48.2579*u.deg, lon=17.0272*u.deg, height=208.0*u.m)
# Set tle
tle = ["ISSd",
"1 25544U 98067A 04160.42390752 .00014992 00000-0 13290-3 0 9491",
"2 25544 51.6329 10.4117 0005395 206.7073 225.7658 15.68815833316945"]
# Set satellite
satellite = twoline2rv(tle[1], tle[2], wgs84)
p = np.zeros(len(dts)*3).reshape(3, len(dts)).astype("float32")
v = np.zeros(len(dts)*3).reshape(3, len(dts)).astype("float32")
for i, tm in enumerate(tms):
p[:, i], v[:, i] = satellite.propagate(tm.datetime.year, tm.datetime.month, tm.datetime.day, tm.datetime.hour, tm.datetime.minute, tm.datetime.second)
# Nutation and precession
epsa, dpsi, deps = nutation_components2000B(t.jd)
rp = precession_matrix_Capitaine(t, Time("J2000"))
rn = nutation_matrix(t)
re = rotation_matrix(-dpsi*np.cos(epsa), "z")
r = matrix_product(rp, rn, re)
pj2000 = np.dot(r, np.array(p))
vj2000 = np.dot(r, np.array(v))
satpos = CartesianRepresentation(x=pj2000[0], y=pj2000[1], z=pj2000[2], unit=u.km)
obspos, obsvel = loc.get_gcrs_posvel(obstime=t)
dr = satpos-obspos
dist, dec, ra = cartesian_to_spherical(dr.x, dr.y, dr.z)
pos = SkyCoord(ra=ra, dec=dec, frame="gcrs")
print('ISS SkyCoord (GCRS) - sgp4 + astropy')
print('{:19s}\t{:7s}\t{:7s}\t{:7s}'.format('Time', 'RA/°', 'DEC/°', 'r/km'))
for tm, p, km in zip(tms, pos, dist):
print('{}\t{:.2f}\t{:.2f}\t{:.3f}'.format(tm.strftime('%Y-%m-%dT%H:%M:%S'),
p.ra.value,
p.dec.value,
km.value))
# Load solary system ephemeris
solar_system_ephemeris.set("de432s")
# Get position of Venus
pvenus = get_body("venus", t, loc)
avenus = np.arcsin(6051.8*u.km/pvenus.distance.to(u.km))
# Get position of sun
psun = get_body("sun", t, loc)
asun = np.arcsin(u.solRad/psun.distance.to(u.km))
# Set up wcs
w = wcs.WCS(naxis=2)
w.wcs.crval = [psun.ra.degree, psun.dec.degree]
w.wcs.crpix = np.array([0.0, 0.0])
w.wcs.cd = np.array([[1.0, 0.0], [0.0, 1.0]])
w.wcs.ctype = ["RA---TAN", "DEC--TAN"]
# Offset position of venus
rxv, ryv = w.wcs_world2pix(pvenus.ra.degree, pvenus.dec.degree, 1)
# Offset position of ISS
rxs, rys = w.wcs_world2pix(pos.ra.degree, pos.dec.degree, 1)
# Create figure
fig, ax = plt.subplots()
ax.set_aspect(1.0)
# Plot Sun
ax.add_artist(Circle((0.0, 0.0), asun.to(u.deg).value, color="y"))
# Plot Venus
ax.add_artist(Circle((rxv, ryv), avenus.to(u.deg).value, color="k"))
# Plot ISS
ax.plot(rxs, rys, "gray")
ax.set_xlim(0.3, -0.3)
ax.set_ylim(-0.3, 0.3)
ax.grid(alpha=0.5)
ax.set_xlabel("R.A. offset (deg)")
ax.set_ylabel("Decl. offset (deg)")
plt.tight_layout()
plt.savefig("transit_astropy.png", bbox_inches="tight")
Raw
iss_transit_astropy2.py
#!/usr/bin/env python3
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Circle
import astropy.units as u
from astropy.time import Time
from astropy.coordinates import get_body, solar_system_ephemeris
from astropy.coordinates import EarthLocation, GCRS, FK5
from astropy.coordinates import SkyCoord, PrecessedGeocentric, CartesianRepresentation
from astropy.wcs import wcs
from astropy.coordinates.earth_orientation import precession_matrix_Capitaine, nutation_components2000B
from astropy.coordinates.matrix_utilities import rotation_matrix, matrix_product, matrix_transpose
from astropy.coordinates import cartesian_to_spherical
from sgp4.earth_gravity import wgs84
from sgp4.io import twoline2rv
# From astropy, but with units
def nutation_matrix(epoch):
"""
Nutation matrix generated from nutation components.
Matrix converts from mean coordinate to true coordinate as
r_true = M * r_mean
"""
# TODO: implement higher precision 2006/2000A model if requested/needed
epsa, dpsi, deps = nutation_components2000B(epoch.jd) # all in radians
epsa *= u.rad
dpsi *= u.rad
deps *= u.rad
return matrix_product(rotation_matrix(-(epsa + deps), 'x'),
rotation_matrix(-dpsi, 'z'),
rotation_matrix(epsa, 'x'))
if __name__ == "__main__":
# Set time
t = Time("2004-06-08T10:09:17", scale="utc", format="isot")
dts = np.arange(5)*u.s
# Generate timestamps
tms = list([t+dt for dt in dts])
# Set location
loc = EarthLocation(lat=48.2579*u.deg, lon=17.0272*u.deg, height=208.0*u.m)
# Set tle
tle = ["ISSd",
"1 25544U 98067A 04160.42390752 .00014992 00000-0 13290-3 0 9491",
"2 25544 51.6329 10.4117 0005395 206.7073 225.7658 15.68815833316945"]
# Set satellite
satellite = twoline2rv(tle[1], tle[2], wgs84)
ps = np.zeros((len(dts), 3), dtype=np.float32)
vs = np.zeros((len(dts), 3), dtype=np.float32)
for i, tm in enumerate(tms):
# TEME position [x,y,z] & velocity [dx,dy,dz] in km and km/s
# TODO: Replace this method. It has a horrible interface (i.e. I'd like to pass the datetime object at once)!
ps[i],vs[i] = satellite.propagate(tm.datetime.year,
tm.datetime.month,
tm.datetime.day,
tm.datetime.hour,
tm.datetime.minute,
tm.datetime.second)
###################################################
# Use skyfield to convert from TEME to ITRF to GCRF
###################################################
from skyfield.constants import AU_KM, DAY_S
from skyfield.sgp4lib import TEME_to_ITRF, ITRF_to_GCRS2
from skyfield.api import load
skyfield_ts = load.timescale()
satposs = []
rGCRS = np.zeros((len(dts),3), dtype=np.float32)
vGCRS = np.zeros((len(dts),3), dtype=np.float32)
for i, (t, p, v) in enumerate(zip(tms, ps, vs)):
skyfield_t = skyfield_ts.from_astropy(t)
rITRF, vITRF = TEME_to_ITRF(skyfield_t.ut1, p, v)
rGCRS[i], vGCRS[i] = ITRF_to_GCRS2(skyfield_t, rITRF, vITRF)
###################################################
# End of ITRF->GCRF conversion code
###################################################
satpos = CartesianRepresentation(x=rGCRS.T[0], y=rGCRS.T[1], z=rGCRS.T[2], unit=u.km)
obspos, obsvel = loc.get_gcrs_posvel(obstime=t)
dr = satpos-obspos
dist, dec, ra = cartesian_to_spherical(dr.x, dr.y, dr.z)
pos = SkyCoord(ra=ra, dec=dec, frame="gcrs")
print('ISS SkyCoord (GCRS) - sgp4 + skyfield.sgp4lib (for frame conversions)')
print('{:19s}\t{:7s}\t{:7s}\t{:7s}'.format('Time', 'RA/°', 'DEC/°', 'r/km'))
for tm, p, km in zip(tms, pos, dist):
print('{}\t{:.2f}\t{:.2f}\t{:.3f}'.format(tm.strftime('%Y-%m-%dT%H:%M:%S'),
p.ra.value,
p.dec.value,
km.value))
# Load solary system ephemeris
solar_system_ephemeris.set("de432s")
# Get position of Venus
pvenus = get_body("venus", t, loc)
avenus = np.arcsin(6051.8*u.km/pvenus.distance.to(u.km))
# Get position of sun
psun = get_body("sun", t, loc)
asun = np.arcsin(u.solRad/psun.distance.to(u.km))
# Set up wcs
w = wcs.WCS(naxis=2)
w.wcs.crval = [psun.ra.degree, psun.dec.degree]
w.wcs.crpix = np.array([0.0, 0.0])
w.wcs.cd = np.array([[1.0, 0.0], [0.0, 1.0]])
w.wcs.ctype = ["RA---TAN", "DEC--TAN"]
# Offset position of venus
rxv, ryv = w.wcs_world2pix(pvenus.ra.degree, pvenus.dec.degree, 1)
# Offset position of ISS
rxs, rys = w.wcs_world2pix(pos.ra.degree, pos.dec.degree, 1)
# Create figure
fig, ax = plt.subplots()
ax.set_aspect(1.0)
# Plot Sun
ax.add_artist(Circle((0.0, 0.0), asun.to(u.deg).value, color="y"))
# Plot Venus
ax.add_artist(Circle((rxv, ryv), avenus.to(u.deg).value, color="k"))
# Plot ISS
ax.plot(rxs, rys, "gray")
ax.set_xlim(0.3, -0.3)
ax.set_ylim(-0.3, 0.3)
ax.grid(alpha=0.5)
ax.set_xlabel("R.A. offset (deg)")
ax.set_ylabel("Decl. offset (deg)")
plt.tight_layout()
plt.savefig("transit_astropy2.png", bbox_inches="tight")
Raw
iss_transit_skyfield.py
#!/usr/bin/env python3
import numpy as np
from skyfield.api import Topos, load
from skyfield.api import EarthSatellite
import matplotlib.pyplot as plt
from matplotlib.patches import Circle
import astropy.units as u
from astropy.wcs import wcs
if __name__ == "__main__":
tle = ["ISSd",
"1 25544U 98067A 04160.42390752 .00014992 00000-0 13290-3 0 9491",
"2 25544 51.6329 10.4117 0005395 206.7073 225.7658 15.68815833316945"]
satellite = EarthSatellite(tle[1], tle[2], tle[0])
ts = load.timescale()
t = ts.utc(2004, 6, 8, 10, 9, 17)
tsat = ts.utc(2004, 6, 8, 10, 9, range(17, 22))
location = Topos(latitude_degrees=48.2579, longitude_degrees=17.0272, elevation_m=208.0)
ra, dec, dist = (satellite-location).at(tsat).radec() #(epoch='date')
# satellite position in ICRF("J2000"), https://rhodesmill.org/skyfield/earth-satellites.html
# print("Satellite ra/dec/dist: {}"ra, dec, dist.km)
print('ISS position in ICRF("J2000") - skyfield')
print('{:19s}\t{:7s}\t{:7s}\t{:7s}'.format('Time', 'RA/°', 'DEC/°', 'r/km'))
for t, ra_deg, dec_deg, km in zip(tsat, ra.to(u.deg).value, dec.to(u.deg).value, dist.km):
print('{}\t{:.2f}\t{:.2f}\t{:.3f}'.format(t.utc_strftime('%Y-%m-%dT%H:%M:%S'),
ra_deg,
dec_deg,
km))
planets = load('de421.bsp')
sun = planets['sun']
earth = planets['earth']
venus = planets['venus']
psun = earth.at(t).observe(sun).apparent().radec()
asun = np.arcsin(u.solRad.to(u.km)/psun[2].km)*u.rad
pvenus = earth.at(t).observe(venus).apparent().radec()
avenus = np.arcsin(6051.8/pvenus[2].km)*u.rad
# Set up wcs
w = wcs.WCS(naxis=2)
w.wcs.crval = [psun[0].to(u.deg).value, psun[1].to(u.deg).value]
w.wcs.crpix = np.array([0.0, 0.0])
w.wcs.cd = np.array([[1.0, 0.0], [0.0, 1.0]])
w.wcs.ctype = ["RA---TAN", "DEC--TAN"]
# Offset position of venus
rxv, ryv = w.wcs_world2pix(pvenus[0].to(u.deg).value, pvenus[1].to(u.deg).value, 1)
# Offset position of ISS
rxs, rys = w.wcs_world2pix(ra.to(u.deg).value, dec.to(u.deg).value, 1)
# Create figure
fig, ax = plt.subplots()
ax.set_aspect(1.0)
# Plot Sun
ax.add_artist(Circle((0.0, 0.0), asun.to(u.deg).value, color="y"))
# Plot Venus
ax.add_artist(Circle((rxv, ryv), avenus.to(u.deg).value, color="k"))
# Plot ISS
ax.plot(rxs, rys, "gray")
ax.set_xlim(0.3, -0.3)
ax.set_ylim(-0.3, 0.3)
ax.grid(alpha=0.5)
ax.set_xlabel("R.A. offset (deg)")
ax.set_ylabel("Decl. offset (deg)")
plt.tight_layout()
plt.savefig("transit_skyfield.png", bbox_inches="tight")
Raw
requirements.txt
astropy==3.2.3
matplotlib==3.1.2
numpy==1.17.4
sgp4==1.4
skyfield==1.15
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment