Imaging observations of HL Tau

The purpose of this notebook is to show how to use a FITS image as input to create a Source object for imaging observations with METIS. For L-band we will work with an ALMA 233 GHz image of HL Tau, pretending that this is the structure of the object in the mid-infrared. For the N-band, we will use a simulation of an AGN.

import os
import numpy as np
from astropy.io import fits
from astropy import units as u
from astropy.wcs import WCS

from matplotlib import pyplot as plt
from matplotlib.colors import LogNorm
%matplotlib inline

import scopesim as sim
sim.bug_report()

# Edit this path if you have a custom install directory, otherwise comment it out. [For ReadTheDocs only]
sim.rc.__config__["!SIM.file.local_packages_path"] = "../../../"

To simulate observations with METIS, the instrument packages for METIS, the ELT, and Armazones are required. The packages are downloaded from the server and installed into sub-directory inst_pkgs in the current working directory. If you have not got the packages yet, uncomment the following cell.

# sim.download_package(['instruments/METIS', 'telescopes/ELT', 'locations/Armazones'])

This notebook uses two example files that can be downloaded from our server. If you already have these files, make sure that they are located in the current working directory. You can then skip the next cell.

sim.download_example_data(["HL_Tau_prep_for_Scopesim.fits", "AGN_uc0890_image_l12_i090_p000.fits"])

Creating a source object from an image

The image that we use for the L-band is HL_Tau_prep_for_Scopesim.fits. This has been derived from HLTau_B6cont_mscale_ap.image.fits, which was retrieved from the ESO Science archive. The main step in preparing this file was a rescaling of the image values such that they directly represent flux in Jy per pixel (the flux scale was arbitrarily set such that the average over the first ring corresponds to 0.01 Jy/arcsec2). The image has a pixel scale of 5 mas, which was retained.

input_hdul = fits.open("HL_Tau_prep_for_Scopesim.fits")

The header contains the information necessary for Scopesim, a WCS and the BUNIT keyword that gives the units of the pixel values. Scopesim mainly needs the CDELT keywords, which contain the pixel scale of the input image. The reference coordinates CRVAL have been set to zero as other values may confuse Scopesim at this stage.

input_hdul[1].header

There are several ways to instantiate a Source object in ScopeSim. In the present case, we provide an image (given as a FITS HDU). The header keyword BUNIT gives the units of the pixel values.

src = sim.Source(image_hdu=input_hdul[1])

L-band observations of the source

The source as defined above could be observed in any waveband – we start in the L-band. The configuration is defined by instantiating a UserCommands object:

cmd_l = sim.UserCommands(use_instrument="METIS", set_modes=["img_lm"])

The default filter is L’:

cmd_l["!OBS.filter_name"]

'Lp'

A different filter could be selected by setting this parameter to a different value. Parameters can also be set by giving a properties dictionary, e.g:

cmd_l = sim.UserCommands(use_instrument="METIS", set_modes=["img_lm"],
                         properties={"!OBS.filter_name": "Mp", "!OBS.exptime": 100.})

We shall stick to the L’ filter here, and set the exposure time directly in the readout() commands later.

We will now build the OpticalTrain (which includes atmosphere, telescope and instrument) and observe the source.

metis_l = sim.OpticalTrain(cmd_l)

metis_l.observe(src, update=True)

The observe command creates the ideal image that arrives at the detector focal plane and stores it as

metis_l.image_planes[0]

This is an HDU object, the image is accessed by appending .data.

Photon noise and detector noise are created by the readout command, which can be applied several times. In this example we read out with two different exposure times. The exposure time is split automatically into DITs such that the detector is not saturated.

result_10 = metis_l.readout(exptime=10.)[0]
print("----")
result_100 = metis_l.readout(exptime=100.)[0]

readout creates a list of HDU lists (a list of multi-extension fits files). In this example the list contains a single HDU list, which we have already extracted by adding [0] above. The detector image is in the first extension of the HDU list. In the following cells we display (sections of) the detector images, with the input image for comparison, and then save the results to fits files.

plt.figure(figsize=(15, 5))
plt.subplot(131)
plt.imshow(input_hdul[1].data[520:1160, 483:1123], origin='lower', vmin=0, vmax=5e-7,
          extent=(-1.6, 1.6, -1.6, 1.6))
plt.xlabel("[arcsec]")
plt.ylabel("[arcsec]")
plt.title("Input image")
plt.subplot(132)
plt.imshow(result_10[1].data[768:1352, 732:1316], origin='lower', extent=(-1.6, 1.6, -1.6, 1.6))
plt.xlabel("[arcsec]")
plt.ylabel("[arcsec]")
plt.title(r"Scopesim simulation, $T_\mathrm{exp} = 10\,\mathrm{s}$")
plt.subplot(133)
plt.imshow(result_100[1].data[768:1352, 732:1316], origin='lower', extent=(-1.6, 1.6, -1.6, 1.6))
plt.xlabel("[arcsec]")
plt.ylabel("[arcsec]")
plt.title(r"Scopesim simulation, $T_\mathrm{exp} = 100\,\mathrm{s}$");

result_10.writeto("hl_tau_img_l_10s.fits", overwrite=True)
result_100.writeto("hl_tau_img_l_100s.fits", overwrite=True)

N-band observation of the source

First we load the image, here an image from a simulation of an AGN with 1% Eddington ratio, observed edge-on at 12 micron.

input_hdul = fits.open("AGN_uc0890_image_l12_i090_p000.fits")

Let’s have a look at the image’s header to see if it is complete for ScopeSim:

input_hdul[0].header

It is complete, but the pixel scale is very small, since the image was created for VLTI observations. We modify the pixel scale and also adjust the total flux to represent a real object, in this case a powerful AGN with a total extent of about 1” and a total flux of about 1 Jy.

input_hdul[0].header["CDELT1"]=1.39e-6
input_hdul[0].header["CDELT2"]=1.39e-6 ## increase pixel scaling by factor of 10
input_hdul[0].data*=50 ## sets total flux to ca. 1 Jy

Let us confirm that the total flux is now about 1 Jy:

np.sum(input_hdul[0].data)

We now set up a new source object. Unfortunately, the value ‘JY/PIXEL’ of the header keyword BUNIT is not compliant with the FITS standard and is not understood by Scopesim. We could simply change it to the compliant form ‘Jy’. Alternatively, the flux parameter can be used to stipulate that an image value of 1 corresponds to 1 Jy. Implicitly, this assigns a spectrum that has a constant value of \(f_{\nu}\) to each pixel.

src = sim.Source(image_hdu=input_hdul[0], flux=1*u.Jy)

We set up the instrument for N-band imaging, the total requested exposure time is set to 3600 s.

cmd_n = sim.UserCommands(use_instrument='METIS', set_modes=['img_n'],
                        properties={"!OBS.exptime": 3600})

The default filter is N2. As before it could be changed by modifying the configuration keyword

cmd_n['!OBS.filter_name']

The instrument itself (the OpticalTrain), including atmosphere and telescope, is built by

metis_n = sim.OpticalTrain(cmd_n)

The list of “effects” in the OpticalTrain can be accessed by metis_n.effects (the method show_in_notebook adds some eye candy for display in this jupyter notebook). The effects that are actually applied are those with included=True, the others are ignored. Note that there are effects that serve the same purpose albeit in different ways. For example, both armazones_atmo_skycalc_ter_curve and armazones_atmo_default_ter_curve provide atmospheric absorption and emission. The first effect downloads the data from ESO’s skycalc server; the second effect uses default data included in the instrument package.

metis_n.effects.show_in_notebook()

We can now “observe” the source. As before, this command creates the ideal image just in front of the detector.

metis_n.observe(src, update=True)

The N band configuration of METIS includes a detector effect ChopNodCombiner. By default, this effect is switched off so that the readout corresponds to a staring observation, which should be sufficient for most science applications. When the effect is switched on it creates four detector images with chopping and nodding offsets and combines them into a chop-nod image with positive and negative beam images. This corresponds to the primary data product of the METIS data reduction pipeline. The default pattern is perpendicular chopping and nodding with

metis_n['chop_nod'].include = True
print("Chopping:", metis_n.cmds['!OBS.chop_offsets'])
print("Nodding: ", metis_n.cmds['!OBS.nod_offsets'])

The exposure time given to the readout() method is the integration time on one position. With exptime=3600, we therefore get a total integration time of 14400 s:

outhdul_perpendicular = metis_n.readout(exptime=3600)[0]

ndit = metis_n.cmds['!OBS.ndit']
dit = metis_n.cmds['!OBS.dit']
outhdul_perpendicular[0].header['DIT'] = dit
outhdul_perpendicular[0].header['NDIT'] = ndit * 4
outhdul_perpendicular[0].header['EXPTIME'] = 4 * dit * ndit
outhdul_perpendicular.writeto('agn_er01_90deg_n_img_perpendicular.fits', overwrite=True)

fig = plt.figure(figsize=(10, 8))
plt.imshow(outhdul_perpendicular[1].data[750:1750, 750:1750], origin='lower', extent=(750, 1750, 750, 1750))
plt.colorbar();

The following cell combines the four beams into a single image. The chop/nod offsets were 3 arcsec, the N-band pixel scale is 6.79 mas.

image = outhdul_perpendicular[1].data
shift = int(3 / 0.00679)
xmin, xmax, ymin, ymax = 823, 1232, 853, 1253
beam_1 = image[ymin:ymax, xmin:xmax]
beam_2 = image[ymin:ymax, xmin+shift:xmax+shift]
beam_3 = image[ymin+shift:ymax+shift, xmin:xmax]
beam_4 = image[ymin+shift:ymax+shift, xmin+shift:xmax+shift]
final_img = beam_1 - beam_2 - beam_3 + beam_4

plt.figure(figsize=(10, 8))
hw_arcsec = (xmax - xmin)/2 * 0.00679
plt.imshow(final_img, origin='lower', extent=(-hw_arcsec, hw_arcsec, -hw_arcsec, hw_arcsec))
plt.xlabel("[arcsec]")
plt.ylabel("[arcsec]")
plt.colorbar();

final_wcs = WCS(naxis=2)
final_wcs.wcs.ctype = ['LINEAR', 'LINEAR']
final_wcs.wcs.crpix = [(final_img.shape[1] - 1)/2, (final_img.shape[0] - 1)/2]
final_wcs.wcs.crval = [0., 0.]
final_wcs.wcs.cdelt = [0.00679, 0.00679] ## this is the METIS detector plate scale in arcseconds
final_wcs.wcs.cunit = ['arcsec', 'arcsec']
final_hdr = final_wcs.to_header()
final_hdr['EXPTIME'] = 4 * dit * ndit
fits.writeto("agn_er01_90deg_n_img_reconstructed.fits", data=final_img, header=final_hdr, overwrite=True)