Setting the correct position for SOHO in a LASCO C3 Map

How to get the correct location of SOHO using JPL HORIZONS and update the header.

This requires the intallation of the astroquery package and an internet connection. astroquery can be installed ontop of the existing sunpy conda environemnt: conda install -c astropy astroquery

# sphinx_gallery_thumbnail_number = 2
import numpy as np
import matplotlib.pyplot as plt

import sunpy.map
from sunpy.coordinates.ephemeris import get_horizons_coord, get_body_heliographic_stonyhurst
from sunpy.net import helioviewer
hv = helioviewer.HelioviewerClient()

Let’s download a SOHO/LASCO C3 image from Helioviewer which provided pre-processed images and load it into a Map.

f = hv.download_jp2('2000/02/27 07:42', observatory='SOHO', instrument='LASCO', detector='C3')
lasco = sunpy.map.Map(f)

Out:

Files Downloaded:   0%|          | 0/1 [00:00<?, ?file/s]

2000_02_27__07_42_05_810__SOHO_LASCO_C3_white-light.jp2:   0%|          | 0.00/528k [00:00<?, ?B/s]

                                                                                                   
Files Downloaded: 100%|##########| 1/1 [00:00<00:00,  6.13file/s]
Files Downloaded: 100%|##########| 1/1 [00:00<00:00,  6.11file/s]
/home/docs/checkouts/readthedocs.org/user_builds/sunpy/conda/stable/lib/python3.7/asyncio/base_events.py:623: ResourceWarning: unclosed event loop <_UnixSelectorEventLoop running=False closed=False debug=False>
  source=self)
/home/docs/checkouts/readthedocs.org/user_builds/sunpy/conda/stable/lib/python3.7/importlib/_bootstrap_external.py:525: ResourceWarning: unclosed <socket.socket fd=16, family=AddressFamily.AF_INET, type=SocketKind.SOCK_STREAM, proto=6, laddr=('172.17.0.3', 53998), raddr=('129.164.137.18', 443)>
  code = marshal.loads(data)
/home/docs/checkouts/readthedocs.org/user_builds/sunpy/conda/stable/lib/python3.7/asyncio/selector_events.py:655: ResourceWarning: unclosed transport <_SelectorSocketTransport fd=16>
  source=self)
/home/docs/checkouts/readthedocs.org/user_builds/sunpy/conda/stable/lib/python3.7/site-packages/glymur/jp2k.py:1171: DeprecationWarning: Use array-style slicing instead.
  warnings.warn("Use array-style slicing instead.", DeprecationWarning)

A user warning let’s you know that there is missing metadata for the observer location. SunPy goes ahead and assumes that the observer is at Earth.

print(lasco.observer_coordinate)

Out:

/home/docs/checkouts/readthedocs.org/user_builds/sunpy/conda/stable/lib/python3.7/site-packages/sunpy-1.0.3-py3.7-linux-x86_64.egg/sunpy/map/mapbase.py:601: SunpyUserWarning: Missing metadata for observer: assuming Earth-based observer.The following sets of keys were checked:
('hgln_obs', 'hglt_obs', 'dsun_obs')
('crln_obs', 'crlt_obs', 'dsun_obs')
  warnings.warn(warning_message, SunpyUserWarning)
<SkyCoord (HeliographicStonyhurst: obstime=2000-02-27T07:42:05.810): (lon, lat, radius) in (deg, deg, km)
    (0., -7.18539124, 1.48139588e+08)>

To check that this worked let’s get the location of Mercury in this exposure and show that it is in the correct location.

mercury_wrong = get_body_heliographic_stonyhurst('mercury', lasco.date, observer=lasco.observer_coordinate)
mercury_hpc_wrong = mercury_wrong.transform_to(lasco.coordinate_frame)
print(mercury_hpc_wrong)

Out:

/home/docs/checkouts/readthedocs.org/user_builds/sunpy/conda/stable/lib/python3.7/site-packages/sunpy-1.0.3-py3.7-linux-x86_64.egg/sunpy/map/mapbase.py:601: SunpyUserWarning: Missing metadata for observer: assuming Earth-based observer.The following sets of keys were checked:
('hgln_obs', 'hglt_obs', 'dsun_obs')
('crln_obs', 'crlt_obs', 'dsun_obs')
  warnings.warn(warning_message, SunpyUserWarning)
INFO: Apparent body location accounts for 330.72 seconds of light travel time [sunpy.coordinates.ephemeris]
/home/docs/checkouts/readthedocs.org/user_builds/sunpy/conda/stable/lib/python3.7/site-packages/sunpy-1.0.3-py3.7-linux-x86_64.egg/sunpy/map/mapbase.py:601: SunpyUserWarning: Missing metadata for observer: assuming Earth-based observer.The following sets of keys were checked:
('hgln_obs', 'hglt_obs', 'dsun_obs')
('crln_obs', 'crlt_obs', 'dsun_obs')
  warnings.warn(warning_message, SunpyUserWarning)
<Helioprojective Coordinate (obstime=2000-02-27T07:42:05.810, rsun=695700000.0 m, observer=<HeliographicStonyhurst Coordinate (obstime=2000-02-27T07:42:05.810): (lon, lat, radius) in (deg, deg, km)
    (0., -7.18539124, 1.48139588e+08)>): (Tx, Ty, distance) in (arcsec, arcsec, km)
    (-23270.61112779, 13370.26466196, 99146796.14617452)>

Let’s plot how this looks with the incorrect observer information.

ax = plt.subplot(projection=lasco)

# Let's tweak the axis to show in degrees instead of arcsec
lon, lat = ax.coords
lon.set_major_formatter('d.dd')
lat.set_major_formatter('d.dd')
ax.plot_coord(mercury_hpc_wrong, 's', color='white', fillstyle='none', markersize=12, label='Mercury')
lasco.plot()
plt.show()
../../../_images/sphx_glr_getting_lasco_observer_location_001.png

Out:

/home/docs/checkouts/readthedocs.org/user_builds/sunpy/conda/stable/lib/python3.7/site-packages/sunpy-1.0.3-py3.7-linux-x86_64.egg/sunpy/map/mapbase.py:601: SunpyUserWarning: Missing metadata for observer: assuming Earth-based observer.The following sets of keys were checked:
('hgln_obs', 'hglt_obs', 'dsun_obs')
('crln_obs', 'crlt_obs', 'dsun_obs')
  warnings.warn(warning_message, SunpyUserWarning)
/home/docs/checkouts/readthedocs.org/user_builds/sunpy/conda/stable/lib/python3.7/site-packages/sunpy-1.0.3-py3.7-linux-x86_64.egg/sunpy/map/mapbase.py:601: SunpyUserWarning: Missing metadata for observer: assuming Earth-based observer.The following sets of keys were checked:
('hgln_obs', 'hglt_obs', 'dsun_obs')
('crln_obs', 'crlt_obs', 'dsun_obs')
  warnings.warn(warning_message, SunpyUserWarning)

SOHO is actually a halo orbit around the Sun–Earth L1 point, about 1 million km away from the Earth. The following functions queries JPL HORIZONS which includes positions of major spacecraft. This function requires an internet connection to fetch the ephemeris data.

soho = get_horizons_coord('SOHO', lasco.date)

Out:

/home/docs/checkouts/readthedocs.org/user_builds/sunpy/conda/stable/lib/python3.7/site-packages/astropy/extern/six.py:7: DeprecationWarning: the imp module is deprecated in favour of importlib; see the module's documentation for alternative uses
  import imp
/home/docs/checkouts/readthedocs.org/user_builds/sunpy/conda/stable/lib/python3.7/site-packages/astropy/extern/six.py:15: AstropyDeprecationWarning: astropy.extern.six will be removed in 4.0, use the six module directly if it is still needed
  AstropyDeprecationWarning)
/home/docs/checkouts/readthedocs.org/user_builds/sunpy/conda/stable/lib/python3.7/site-packages/astropy/extern/bundled/six.py:66: ResourceWarning: unclosed file <_io.TextIOWrapper name='/home/docs/checkouts/readthedocs.org/user_builds/sunpy/conda/stable/lib/python3.7/site-packages/astropy/extern/bundled/six.py' mode='r' encoding='utf-8'>
  except OverflowError:
INFO: Obtained JPL HORIZONS location for SOHO (spacecraft) (-21) [sunpy.coordinates.ephemeris]

For fun, let’s see how far away from the Earth SOHO is by converting to an Earth-centered coordinate system (GCRS).

print(soho.transform_to('gcrs').distance.to('km'))

Out:

1352066.1922767868 km

Let’s fix the header.

lasco.meta['HGLN_OBS'] = soho.lon.to('deg').value
lasco.meta['HGLT_OBS'] = soho.lat.to('deg').value
lasco.meta['DSUN_OBS'] = soho.radius.to('m').value

Let’s get the right position now.

mercury = get_body_heliographic_stonyhurst('mercury', lasco.date, observer=lasco.observer_coordinate)
mercury_hpc = mercury.transform_to(lasco.coordinate_frame)

Out:

INFO: Apparent body location accounts for 326.33 seconds of light travel time [sunpy.coordinates.ephemeris]

The difference between the incorrect position and the right one is:

r = np.sqrt((mercury_hpc.Tx - mercury_hpc_wrong.Tx) ** 2 + (mercury_hpc.Ty - mercury_hpc_wrong.Ty) ** 2)
print(r)

Out:

136.73644924312012 arcsec

Let’s plot the results.

ax = plt.subplot(projection=lasco)

# Let's tweak the axis to show in degrees instead of arcsec
lon, lat = ax.coords
lon.set_major_formatter('d.dd')
lat.set_major_formatter('d.dd')
ax.plot_coord(mercury_hpc, 's', color='white', fillstyle='none', markersize=12, label='Mercury')
lasco.plot()
plt.show()
../../../_images/sphx_glr_getting_lasco_observer_location_002.png

Total running time of the script: ( 0 minutes 1.967 seconds)

Gallery generated by Sphinx-Gallery