Getting HST Ephemeris and Related Data for an Observation

Worked Examples for an ACS/SBC Observation

By Chris Clark (ESA/AURA, Space Telescope Science Institute; cclark@stsci.edu), last edited 28th March 2025

This notebook walks through how to get various information about where HST is, and where it is pointing, in relation to Earth and the Sun, for a given observation.

The methods laid out in this notebook do not necessarily represent the “official” way that one would extract these quantities from the full NASA engineering data. However, interacting with the engineering data can be non-trivial, and there is often limited documentation.

Therefore, this notebook aims to provide ways to calculate various ephemeris (and related) quantities in ways that are reliable, whilst also still being relatively straightforward.

Some of the examples laid out in this notebook use information available in the FITS header, and some use the JPL Horizons service or pvlib package. To run this notebook, you will require astropy, astroquery, and pvlib to be installed (along with their dependences).

In this example, we will be using an observation from the ACS/SBC - namely jec0c6aeq_flc.fits. It’s an observation of part of nearby galaxy NGC7793, looking at a bunch of young stars and star-forming regions. This file will be downloaded from MAST, and put in a directory your current working directory (in general, this will automatically be the same directory as the notebook itself, under standard notebook kernel server settings).

Note that the methods demonstrated in this notebook require you to be working with a file associated with an individual expoure (ie, a raw, flt, or flc file). This is because a drizzled file, or other combined data product, will consist of data from different observations, where the telescope, Earth, and Sun will have been in different relative positions.

Table of Contents

1. Earth Limb Angle

2. HST Lattude & Longitude During Observation

3. Solar Time, Solar Altitude, and Solar Cycles

# Required imports all here
import os
import numpy as np
import astropy.io.fits
import astropy.wcs
import astropy.time
import astropy.convolution
import astropy.units as u
import astropy.coordinates
import astroquery.mast
from astroquery.jplhorizons import Horizons
import pvlib

# Report to user where data folder will be palced
data_dir = os.path.join(os.getcwd())
print('Data for this notebook will be downloaded to:')
print(data_dir)

Data for this notebook will be downloaded to:
/home/runner/work/hst_notebooks/hst_notebooks/notebooks/ACS/hst_orbits_ephem

Before we get going, we’ll download the obesrvation we want from MAST, and read in its headers.

# Create folder for holding downloaded files, etc

if not os.path.exists(data_dir):
    os.mkdir(data_dir)

# Download our exmaple observation
flt_filename = 'jec0c6aeq_flt.fits'
flt_path = os.path.join(data_dir, flt_filename)
mast_msg = astroquery.mast.Observations.download_file('mast:HST/product/'+flt_filename,
                                                      local_path=flt_path,
                                                      cache=False)
print(mast_msg[0])

# Read in headers for our observation; we will require info from the primary and image headers
flt_hdr_0 = astropy.io.fits.getheader(flt_path, ext=0)
flt_hdr_1 = astropy.io.fits.getheader(flt_path, ext=1)

Downloading URL https://mast.stsci.edu/api/v0.1/Download/file?uri=mast:HST/product/jec0c6aeq_flt.fits to /home/runner/work/hst_notebooks/hst_notebooks/notebooks/ACS/hst_orbits_ephem/jec0c6aeq_flt.fits ...

 [Done]

COMPLETE

1. Earth Limb Angle

Earth Limb Angle (ELA) is the angle between the locaiton where the telescope is pointing, and the limb of the Earth. As ELA gets smaller, the atmosphere can have an increasing effect on HST observations in certain wavelenghts (especially the IR and UV). So knowing the ELA is worthwhile for understanding backgrounds, etc.

Establishing ELA for a given observation is a somewhat involved process, so we will deal with it first.

To get the ELA, we need to get the jitter file (ie, _jit.fits file for this observation. If this observation was not part of an association, then the jit file should have the same rootname (ie, jec0c6aeq) as the flt file. However, if this oabservation was part of an assocation, then the jit file will have the same rootname as the associaiton. So we will check for both

# First, check if jit file exists with same rootname as flt
jit_filename = flt_filename.replace('_flt.fits', '_jit.fits')
jit_path = os.path.join(data_dir, jit_filename)
mast_msg = astroquery.mast.Observations.download_file('mast:HST/product/'+jit_filename,
                                                      local_path=jit_path,
                                                      cache=False)
if os.path.exists(jit_path) and (mast_msg[0] != 'ERROR'):
    jit_data = astropy.io.fits.getdata(jit_path)

# Otherwise, get the jit file by finding the association of this observation
elif 'HTTPError: 404 Client Error: Not Found for url' in mast_msg[1]:
    asn_id = flt_hdr_0['ASN_ID'].lower()
    jit_filename = asn_id+'_jit.fits'
    jit_path = os.path.join(data_dir, jit_filename)
    mast_msg = astroquery.mast.Observations.download_file('mast:HST/product/'+jit_filename,
                                                          local_path=jit_path,
                                                          cache=False)
    jit_data = astropy.io.fits.getdata(jit_path)
print(mast_msg[0])

Downloading URL https://mast.stsci.edu/api/v0.1/Download/file?uri=mast:HST/product/jec0c6010_jit.fits to /home/runner/work/hst_notebooks/hst_notebooks/notebooks/ACS/hst_orbits_ephem/jec0c6010_jit.fits ...

 [Done]

COMPLETE

Now we have the jit file read it, we grab the ELA data from it. This consists of an array of ELA values, as the ELA changes continually throughout the exposure.

# Get ELA data from JIT file
ela_arr = np.array(jit_data['LimbAng'])

# Check that ELA data is good (can be an empty array for certain engineering data, etc)
try:
    ela_min = np.min(ela_arr)
except ValueError:
    print('ELA array is empty')

# Assuming ELA data is good, get some summary statistics
ela_min = np.min(ela_arr)
ela_max = np.max(ela_arr)
ela_mean = np.mean(ela_arr)
print(f'ELA ranges from {ela_min:.1f} deg to {ela_max:.1f} deg, with a mean of {ela_mean:.1f} deg')

ELA ranges from 39.6 deg to 52.9 deg, with a mean of 46.3 deg

2. HST Lattude & Longitude During Observation

If you want to know the position of HST over the Earth during a given observation, you have to account for the fact that HST is continually moving in three dimensions around the Earth, which is continually moving in three dimensions around the Sun (thanks, Copernicus). This is a whole bunch of 3D ephemeris maths we’d rather not do by hand.

Fortunately, we don’t have to! We can use astropy.coordinates and the JPL Horizons service instead.

First, we have to know exactly when the observation was taken, which we can get from the header. We then use that to construct an astropy.time.Time object.

# Get observation date and time fromt the header
date_string = flt_hdr_0['DATE-OBS']
time_string = flt_hdr_0['TIME-OBS']

# Convert into an astropy Time object
datetime_obs_start = astropy.time.Time(date_string+'T'+time_string, format='isot', scale='utc')

Now we use the datetime, and the exposure time, to find the mid-point time of the exposure, and construct an array sampling the whole duration of the observation. We also compute the Modified Julian Date (MJD) of the mid-point of the exposure, as this will be a useful number later

# Use duration of exposure to find time in middle of exposure
exp_time = flt_hdr_0['EXPTIME'] * u.second
datetime_mid_delta = exp_time * 0.5
datetime_obs_mid = datetime_obs_start + datetime_mid_delta

# Make array densely sampling time throughot the observation (1000 points is a good arbitrary dense sampling)
datetime_obs_arr = np.linspace(datetime_obs_start, (datetime_obs_start + exp_time), num=1000)

# Record date and time info
mjd_obs_mid = astropy.time.Time(datetime_obs_mid).mjd

# Report summary of observation time
print(f'Exposure started at {datetime_obs_arr[0]}, and ended at {datetime_obs_arr[-1]}')
print(f'The mid-point time of the exposure was at MJD of {mjd_obs_mid:.5f}')

Exposure started at 2022-04-25T03:18:51.000, and ended at 2022-04-25T03:22:41.000
The mid-point time of the exposure was at MJD of 59694.13942

With all this date & time info ready, we can now use the JPL Horizons service to work out exactly where HST was during this exposure.

We can query the service via astroquery.jplhorizons.Horizons. The Horizons service uses Navigation and Ancillary Information Facility (NAIF) ID codes to identify objects and coordinates. HST has a dedicated NAIF ID code: namely, -48. We set the location of the query to 500@399; this is combination of NAIF ID codes that corresponds to the location of Earth (the @339 part), described in a Geocentric frame (the 500 part). The Horizons and astropy.coordinates systems are smart enough that nothing critically depends on choosing the ‘correct’ location here - but as HST orbits the Earth, this reference frame will make all the numbers involved a little more easily interpretable.

# Query JPL Horizons, getting values for middle of observation
jplhz_hst_id = -48
jplhz_obj_obs_mid = Horizons(id=jplhz_hst_id, location='500@399', epochs=datetime_obs_mid.jd)

Now we extract the 3D position of HST in relation to Earth, in the International Terrestrial Reference System (ITRS), and use this to construct an astropy.coordinates object for Earth’s location.

# Get HST vectors from Horizons object
jplhz_vector = jplhz_obj_obs_mid.vectors()

# Extract position in ITRS (International Terrestrial Reference System)
hst_geo_x = jplhz_vector['x'][0] * u.AU.to(u.m)
hst_geo_y = jplhz_vector['y'][0] * u.AU.to(u.m)
hst_geo_z = jplhz_vector['z'][0] * u.AU.to(u.m)

# Create an EarthLocation object for the HST
hst_location = astropy.coordinates.EarthLocation.from_geocentric(hst_geo_x, hst_geo_y, hst_geo_z, unit=u.m)

In order to use this to work out the terrestrial latitude and longitude that HST was above during observations, we need to convert to geodetic coordinates, and use that to build a astropy.coordinates.EarthLocation object.

# Compute HST location in map coordiantes
hst_location = hst_location.to_geodetic()
hst_location = astropy.coordinates.EarthLocation.from_geodetic(lat=hst_location.lat,
                                                               lon=hst_location.lon,
                                                               height=hst_location.height)

# Report location to user
hst_location_lat = hst_location.lat.value
hst_location_lon = hst_location.lon.value
hst_location_height = hst_location.height.to(u.km)
print(f'During exposure mid-point, HST was above lat = {hst_location_lat:.2f}, lon = {hst_location_lon:.2f}, at {hst_location_height:.0f}')

During exposure mid-point, HST was above lat = -46.67, lon = 46.65, at 547 km

As an example use-case of this, we can compute rough distance from the South Atlantic Anomaly (SAA). Note that the SAA is not symmetric, and that any proper SAA calculations should take account of its contours. And the central point of the SAA changes. So this for illustrative purposes only.

# Give approxiamte coords of centre of SAA
saa_lat, saa_lon = -25.603*u.deg, -40.287*u.deg # Approximate coords, but good enough
saa_location = astropy.coordinates.EarthLocation.from_geodetic(lat=saa_lat, lon=saa_lon)

# Compute separation between SAA and HST, in degrees, and report to user
saa_itrs = astropy.coordinates.SkyCoord(saa_location.get_itrs())
hst_itrs = astropy.coordinates.SkyCoord(hst_location.get_itrs())
saa_distance = saa_itrs.separation(hst_itrs)
print(f'HST was {saa_distance:.1f} from the middle of the SAA')

HST was 69.8 deg from the middle of the SAA

3. Solar Time, Solar Altitude, and Solar Cycles

Now we can do a bunch of Solar maths regarding where, and when, HST peformed this exposure.

To start off, let’s calculate Solar time. This is the time system where the Sun is overhead at 12 noon at whatever your location happens to be. This can be relevant in filters where the atmosphere gets excited by Solar radiation during the course of the day, resulting in increased airglow at dusk vs dawn (even though comparable dust and dawn times may have the same Solar altitude).

We get the Sun’s location using the package pvlib, which is designed for Solar maths. Then we compute the ‘equation of time’ for HST during the exposure. The ‘equation of time’ is the difference between apparent Solar time and mean Solar time for a given time & place.

# First, find the Sun's location 
sun_location = pvlib.location.Location(hst_location.lat.value, hst_location.lon.value).get_solarposition(datetime_obs_mid.datetime)

# Find the equation of time for HST, and thence the longitude correction required to get actual Solar time from HST's position
hst_eot = (sun_location.equation_of_time.array[0] * u.minute).to(u.hour)
hst_lon_corr = hst_location.lon / 15.0 * u.hour/u.deg

# Use the equation fof time to calculate solar time at HST's location, as an astropy.time.Time object
datetime_obs_mid_solar = datetime_obs_mid + hst_eot + hst_lon_corr

# Get out the solar time as an actual number, hours
hst_solar_time_hr = datetime_obs_mid_solar.datetime.hour + (datetime_obs_mid_solar.datetime.minute / 60.0) + (datetime_obs_mid_solar.datetime.second / 3600.0)
print(f'Solar time for HST during exposure was {hst_solar_time_hr:05.2f} hours (in 24-hour notation)')

Solar time for HST during exposure was 06.49 hours (in 24-hour notation)

Now, we’ll calculate the altitude and azimuth of the Sun, from the position of HST. Turns out that astropy.coordinates has a built-in function for this, assuming you construct the right coordinate frame first.

Note that the altitude and azimuth will be for the plane of HST tangent to the surface of the Earth, at the location of HST. Therefore, the Earth limb will be at a negative altitude, due to the height of HST’s orbit above the surface.

# Use astropy coordinates to construct an alt-az reference frame object at the place and time of the exposure
altaz_frame = astropy.coordinates.AltAz(obstime=datetime_obs_mid, location=hst_location)

# Now use astropy function to compute Solar position in this reference frame
sun_altaz = astropy.coordinates.get_sun(datetime_obs_mid).transform_to(altaz_frame)
print(f'Sun\'s position during exposure was altitude = {sun_altaz.alt:.2f}, azimuth = {sun_altaz.az:.2f}')

Sun's position during exposure was altitude = -4.61 deg, azimuth = 75.68 deg

Now to finish with a couple of easy calculations. Firstly, lets find out how far the observations are from the closest Solar cycle maximum (either forward or backward in time). The atmosphere extends to higher altitudes during Solar maxima, so this can once again be a relevant quantitity for airglow.

We take the dates of past maxima from the International Sunspot Numbers v2.0, as available at SILSO. (Maxima dates are reported as months only in ISNv2.0; here we take the middle of each month as the peak for each, for convenience.)

# First, state the dates of of the Solar maxima  the once before HST launch
solar_cycle_peaks_mjd = ['2024-06-15', '2014-04-15', '2001-11-15', '1989-11-15']

# Convert this date strings into astropy.time.Time objects
solar_cycle_peaks_mjd = np.array([astropy.time.Time(solar_cycle_peak_mjd, format='iso').mjd for solar_cycle_peak_mjd in solar_cycle_peaks_mjd])

# Calculate how many days, in MJD, separate our HST observation from the closest solar maximum (forward or backward in time)
solar_peak_gap = np.min(np.abs(mjd_obs_mid - solar_cycle_peaks_mjd))
print(f'Observation happened {solar_peak_gap:.0f} days away from nearest Solar maximum')

Observation happened 782 days away from nearest Solar maximum

Lastly, let’s convert the ra and dec of our HST observations to ecliptic coordinates. This can be useful if we’re worrdied about the potential for zodiacal light to affect our observations. For this calculation, we’re going to grab the centre coord of the image in our flt file.

# Get ra and dec from header; here we just use the WCS reference coord (you may want to be more precise)
ra = flt_hdr_1['CRVAL1A'] * u.deg
dec = flt_hdr_1['CRVAL2A'] * u.deg

# Convert coordinates to ecliptic frame, and report to user
icrs_coords = astropy.coordinates.SkyCoord(ra=ra, dec=dec, frame='icrs')
ecliptic_coords = icrs_coords.transform_to(astropy.coordinates.GeocentricTrueEcliptic)
print(f'Reference pixel of observation has ecliptic coords of lat = {ecliptic_coords.lat:.2f}, lon = {ecliptic_coords.lon:.2f}')

Reference pixel of observation has ecliptic coords of lat = -29.42 deg, lon = 345.26 deg

If you use astropy or pvlib for published research, please cite the authors. Follow these links for more information about citing astropy or pvlib:

• Citing astropy

• Citing pvlib