On May 9, 2016, Mercury made a 7.5 hour transit across the Sun. In this notebook, we will be making a composite image of Mercury's transit across the Sun using data from both the NASA Solar Dynamics Observatory's Helioseismic and Magnetic Imager (HMI) instrument and the ESA Solar and Heliospheric Observatory's Michelson Doppler Imager (MDI) instrument.
import json, urllib, numpy as np, matplotlib.pylab as plt, requests
import pandas as pd
from astropy.io import fits
from matplotlib.cm import register_cmap, cmap_d
from scipy import ndimage
%matplotlib inline
%config InlineBackend.figure_format='retina'
Now we define a color table exactly like the SAO ds9 heat color table:
ds9heat = {'red': lambda v : np.interp(v, [0, 0.34, 1], [0, 1, 1]),
'green': lambda v : np.interp(v, [0, 1], [0, 1]),
'blue': lambda v : np.interp(v, [0, 0.65, 0.98, 1], [0, 0, 1, 1])}
register_cmap('ds9heat', data=ds9heat)
Then we'll grab the data. Normally we pull continuum intensity data from the JSOC database, but these continuum intensity data were specially created from level-1 filtergram data such that the cadence of the images capture Mercury's ingress and egress. For the purpose of a composite image, using these data is a purely aesthetic choice, and the following methodology (with slightly different thresholds) works for the nominal continuum intensity data as well.
url = 'http://jsoc.stanford.edu/data/events/Mercury_HMI_I/WholeSun/fits/'
files = pd.read_html(url)
list_of_files = files[0]['Name'][2:556]
print("There are",len(list_of_files),"files.")
There are 554 files.
The real difficulty in making a composite image is to capture the entire effect of Mercury as the scattered light around the planet creates a couple-pixel annulus with an intensity comparable to other regions of the Sun. To capture this effect, we follow a six-step procedure:
After selecting all of the images that will go into our composite image (about every 20 minutes looks aesthetically pleasing), we subtract from each image a background image during which Mercury is not transiting the Sun. This differenced image (stored in variable differenced_image
) is especially useful in eliminating the effect of limb darkening. We then set all pixels that are not within the transit path to zero.
Then we apply a threshold to the differenced image to capture the pixels that correspond to Mercury. This requires some experiment: if we choose a threshold that is too high, we don't capture the entire planet. However, if we choose a threshold that is too low, we capture extraneous pixels. Therefore, we'll first apply a threshold that is too high -- to make sure we don't capture extraneous pixels -- and then generate a binary mask where the pixels labeled with a 1 represent the location of Mercury and the pixels labeled with a 0 represent the locations where Mercury does not exist (stored in variable mask
).
Now, we'll compensate for the fact that we chose a low threshold: we'll use morphological binary dilation to grow the radius of Mercury within the mask by a few pixels (stored in variable dilated_mask
).
We can create a masked image by multiplying this dilated binary mask by the original image to get the value of the pixels at the locations of Mercury (stored in variable masked_image
).
background = fits.open("http://jsoc.stanford.edu/data/events/Mercury_HMI_I/WholeSun/fits/20160509_185527_UTC.0554.fits")
background = background[1].data
mask_inverse = np.ones([4096,4096])
masked_image = np.zeros([4096,4096])
composite_dilated_mask = np.zeros([4096,4096])
for i in range(14,len(list_of_files),21):
if (i == 539):
i = 540 # capture the egress
mask = np.zeros([4096,4096])
continuum = fits.open("http://jsoc.stanford.edu/data/events/Mercury_HMI_I/WholeSun/fits/"+list_of_files[i])
continuum = continuum[1].data
differenced_image = background - continuum
differenced_image[0:1300,:] = 0.0 # mask out the area outside the transit path
differenced_image[1500:4096,:] = 0.0
find_mercury = np.where(differenced_image > 14000)
for i in range(len(find_mercury[0])):
mask[find_mercury[0][i],find_mercury[1][i]] = 1.
dilated_mask = ndimage.binary_dilation(mask, iterations=2).astype(mask.dtype)
composite_dilated_mask = composite_dilated_mask + dilated_mask
masked_image = continuum*dilated_mask + masked_image
This is what each of these arrays look like within the transit path:
fig = plt.figure(figsize=(14,2.5))
ax = fig.add_subplot(2,2,1)
ax.get_xaxis().set_ticks([])
ax.get_yaxis().set_ticks([])
ax.set_title('differenced image')
plt.imshow(differenced_image[1100:1700,:],cmap='Greys',vmax=30000,vmin=0,origin='lower')
ax = fig.add_subplot(2,2,2)
ax.get_xaxis().set_ticks([])
ax.get_yaxis().set_ticks([])
ax.set_title('mask')
plt.imshow(mask[1100:1700,:],cmap='Greys',vmax=1.,vmin=0,origin='lower')
ax = fig.add_subplot(2,2,3)
ax.get_xaxis().set_ticks([])
ax.get_yaxis().set_ticks([])
ax.set_title('dilated mask')
plt.imshow(dilated_mask[1100:1700,:],cmap='Greys',vmax=1.,vmin=0,origin='lower')
ax = fig.add_subplot(2,2,4)
ax.get_xaxis().set_ticks([])
ax.get_yaxis().set_ticks([])
ax.set_title('masked image')
plt.imshow(masked_image[1100:1700,:],cmap='ds9heat',vmax=20000,vmin=0,origin='lower')
<matplotlib.image.AxesImage at 0x117b1c6a0>
We can then multiply the inverse of the dilated binary mask -- where the the pixels labeled with a 0 represent the location of Mercury and the pixels labeled with a 1 represent the locations where Mercury does not exist (stored in variable inverse_dilated_mask
) -- with the background image to create a composite image (stored in variable composite
) with holes where Mercury should be.
inverse_dilated_mask = np.logical_not(composite_dilated_mask).astype(composite_dilated_mask.dtype)
composite = background * inverse_dilated_mask
fig = plt.figure(figsize=(14,14))
ax = fig.add_subplot(1,1,1)
ax.get_xaxis().set_ticks([])
ax.get_yaxis().set_ticks([])
plt.imshow(composite,cmap='ds9heat',vmax=62000,vmin=0,origin='lower')
<matplotlib.image.AxesImage at 0x118c72860>
Finally, we can fill in these holes with the data by adding the composite image to the masked image.
final = composite + masked_image
print("The maximum value for this array is ",(np.nanmax(final)))
fig = plt.figure(figsize=(14,14))
ax = fig.add_subplot(1,1,1)
ax.get_xaxis().set_ticks([])
ax.get_yaxis().set_ticks([])
plt.imshow(final,cmap='ds9heat',vmax=np.nanmax(final),vmin=0,origin='lower')
The maximum value for this array is 62702.0
<matplotlib.image.AxesImage at 0x105627320>
Finally, save the figure at its native resolution (4096x4096 pixels):
%%capture
# dpi for retina display
# detect monitor dpi automatically here: http://www.infobyip.com/detectmonitordpi.php
my_dpi=192
fig = plt.figure(figsize=(4096.0/my_dpi, 4096.0/my_dpi), dpi=my_dpi)
ax = fig.add_subplot(1,1,1)
fig.subplots_adjust(left=0,right=1,top=1,bottom=0)
ax.get_xaxis().set_ticks([])
ax.get_yaxis().set_ticks([])
plt.imshow(final,cmap='ds9heat',vmax=np.nanmax(final),vmin=0,origin='lower',interpolation='nearest')
fig.savefig('mercury_composite_HMI.jpg')
We can repeat the same process, using different thresholds, with MDI data. Since MDI is an older instrument, the telescope's front window has physically degraded over time. This and other factors contribute to more scattered light within the instrument. As such, the annulus of scattered light around Mercury in the MDI images is much larger than that in the HMI images.
url = 'http://jsoc.stanford.edu/data/events/Mercury_MDI/WholeSun/fits/'
files = pd.read_html(url)
list_of_files = files[0]['Name'][2:150]
print("There are",len(list_of_files),"files.")
There are 148 files.
Some of the WCS keyworks in the MDI FITS files have non-standard values so we use the astropy package to fix these keywords. Also, there is a gap in the transmission of the MDI data at the beginning of the eclipse (not at the ingress, but after the first image); the irregular cadence (rangeofis
) is used to mitigate that effect in the final image.
background = fits.open("http://jsoc.stanford.edu/data/events/Mercury_MDI/WholeSun/fits/20160509_170724_UTC.0148.fits")
background.verify("fix")
background = background[1].data
mask_inverse = np.ones([1024,1024])
masked_image = np.zeros([1024,1024])
composite_dilated_mask = np.zeros([1024,1024])
rangeofis = [0, 4, 6, 12, 18, 24, 30, 36, 42, 48, 54, 60, 66, 72, 78, 84, 90, 96, 102, 108, 114, 120, 126, 132, 137, 142]
for i in rangeofis:
mask = np.zeros([1024,1024])
continuum = fits.open("http://jsoc.stanford.edu/data/events/Mercury_MDI/WholeSun/fits/"+list_of_files[i])
continuum.verify("fix")
continuum = continuum[1].data
find_skewed_pixels = np.where(continuum < 170)
for i in range(len(find_skewed_pixels[0])):
continuum[find_skewed_pixels[0][i],find_skewed_pixels[1][i]] = 0.
differenced_image = background - continuum
differenced_image[0:400,:] = 0.0 # mask out the area outside the transit path
differenced_image[500:1024,:] = 0.0
find_mercury = np.where(differenced_image > 175)
for i in range(len(find_mercury[0])):
mask[find_mercury[0][i],find_mercury[1][i]] = 1.
dilated_mask = ndimage.binary_dilation(mask, iterations=2).astype(mask.dtype)
composite_dilated_mask = composite_dilated_mask + dilated_mask
masked_image = continuum*dilated_mask + masked_image
And finally, we make and save the image. We don't have data during the time of the ingress, but we do see the egress. It is worth noting that the transit path was different from MDI's vantage point at L1 versus HMI's vantage point in inclined geosynchronous orbit.
inverse_dilated_mask = np.logical_not(composite_dilated_mask).astype(composite_dilated_mask.dtype)
composite = background * inverse_dilated_mask
final = composite + masked_image
fig = plt.figure(figsize=(14,14))
ax = fig.add_subplot(1,1,1)
ax.get_xaxis().set_ticks([])
ax.get_yaxis().set_ticks([])
plt.imshow(final,cmap='ds9heat',vmax=700,vmin=200,origin='lower')
<matplotlib.image.AxesImage at 0x106754978>
%%capture
# dpi for retina display
# detect monitor dpi automatically here: http://www.infobyip.com/detectmonitordpi.php
my_dpi=192
fig = plt.figure(figsize=(1024.0/my_dpi, 1024.0/my_dpi), dpi=my_dpi)
ax = fig.add_subplot(1,1,1)
fig.subplots_adjust(left=0,right=1,top=1,bottom=0)
ax.get_xaxis().set_ticks([])
ax.get_yaxis().set_ticks([])
plt.imshow(final,cmap='ds9heat',vmax=700,vmin=200,origin='lower',interpolation='nearest')
fig.savefig('mercury_composite_MDI.jpg')