Jupiter: CIRS nadir-viewing thermal emission
In this example, we show how archNEMESIS can be used to calculate forward models for a nadir-viewing observation. In particular, we are going to model Jupiter’s thermal emission using instrument specifications relevant for the Composite Infrared Spectrometer (CIRS) on Cassini.
First of all, we show how all the input files can be read at once, and summarise the content of these files. Then, we will show how the forward models can be easily computed in python.
[1]:
import archnemesis as ans
import numpy as np
import matplotlib.pyplot as plt
1. Inspecting the input files
archNEMESIS includes the read_input_files() function, which automatically reads all input files and stores the information in the reference classes. Here, we explore some of the most important characteristics of this particular forward model.
[2]:
runname = 'cirstest'
#Reading the input files
Atmosphere,Measurement,Spectroscopy,Scatter,Stellar,Surface,CIA,Layer,Variables,Retrieval = ans.Files.read_input_files(runname)
INFO :: read_input_files :: Files.py-435 :: Reading atmospheric files...
INFO :: read_input_files :: Files.py-576 :: Reading .apr file
INFO :: read_apr :: Variables_0.py-823 ::
Variables_0 :: read_apr :: varident [0 0 0]. Constructed model "Model0" (id=0)
INFO :: read_apr :: Variables_0.py-826 :: Model0:
|- id : 0
|- parent classes: PreRTModelBase
|- description: In this model, the atmospheric parameters are modelled as
| continuous profiles in which each element of the state vector
| corresponds to the atmospheric profile at each altitude level
|- n_state_vector_entries : 81
|- state_vector_slice : slice(0, 81, None)
|- state_vector_start : 0
|- target : 1
|- Parameters:
| |- full_profile :
| | |- slice : slice(None, None, None)
| | |- unit : PROFILE_TYPE
| | |- description: Every value for each level of the profile
| | |- apriori value : [328.16 305.63 284.53 264.74 246.16 228.71 212.31 196.9 182.41 168.83
156.1 144.21 133.41 124.3 117.86 113.89 110.75 108.25 106.26 104.68
103.4 104.5 107.3 110.1 113.5 117. 120.5 125. 128.5 132.4
136. 139.8 143.4 146.2 148.4 150.3 151.8 153. 154.1 155.1
156. 156.83 157.66 158.49 159.32 160.15 160.98 161.81 162.64 163.47
164.3 165.14 165.98 166.82 167.66 168.5 169.34 170.18 171.02 171.86
172.7 173.53 174.36 175.19 176.02 176.85 177.68 178.51 179.34 180.17
181. 185. 189. 193. 197. 201. 208.13 215.27 222.4 229.53
236.67]
Atmosphere
The profiles in this example were identified as standard for CIRS data analysis.
Temperature, pressure and altitude levels provided by Glenn Orton assuming g = 24 m s\(^{-2}\).
Ammonia assumed to be constant in the troposphere, limited by the saturation vapour pressure between 0.5-0.6 bar. Above this level, its abundance drops according to Fouchet et al. (1999).
Phosphine is the standard NIMS case.
Hydrocarbons (including CH\(_4\)) from Romani note to CIRS team (7 May 1999) with deep CH\(_4\) nudged to GPMS results.
H\(_2\) and He from Galileo entry probe.
\(^{13}\)CH\(_4\) (2) scaled from \(^{12}\)CH\(_4\) (1) assuming terrestrial relative abundances.
\(^{12}\)CH\(_3\)D (3) also scaled from \(^{12}\)CH\(_4\) (1) but using \(^{12}\)CH\(_3\)D/\(^{12}\)CH\(_4\) = 8 \(\times\) 10\(^{-5}\).
[3]:
Atmosphere.plot_Atm()
Atmosphere.plot_Dust()
Measurement
In this example, as we only want to model the spectrum, the information about the measured spectra is not relevant. However, the .spx file also defines the geometry of the observations and the spectral range and resolution to use.
The spectral range in this example goes from 5-1399 cm\(^{-1}\). The FWHM of the instrument function is set to 0, which means that it is accounted for directly in the k-tables. In this simple example, all relevant angles were assumed to be zero.
[4]:
Measurement.summary_info()
INFO :: summary_info :: Measurement_0.py-422 :: Spectral resolution of the measurement is account for in the k-tables
INFO :: summary_info :: Measurement_0.py-432 :: Field-of-view centered at :: ('Latitude', 0.0, '- Longitude', 0.0)
INFO :: summary_info :: Measurement_0.py-433 :: There are (1, 'geometries in the measurement vector')
INFO :: summary_info :: Measurement_0.py-435 ::
INFO :: summary_info :: Measurement_0.py-436 :: GEOMETRY 1
INFO :: summary_info :: Measurement_0.py-437 :: Minimum wavelength/wavenumber :: 2000.0 um/5.0 cm^-1 - Maximum wavelength/wavenumber :: 7.147962830593281 um/1399.0 cm^-1
INFO :: summary_info :: Measurement_0.py-470 :: Nadir-viewing geometry. Latitude :: (np.float64(0.0), ' - Longitude :: ', np.float64(0.0), ' - Emission angle :: ', np.float64(0.0), ' - Solar Zenith Angle :: ', np.float64(0.0), ' - Azimuth angle :: ', np.float64(0.0))
Variables
In this example, we choose a parameterisation that allows us to define the temperature in each altitude level. In particular, we choose the model parameterisation #0, which allows the definition of a continuous vertical profile.
[5]:
print('VARIDENT')
print(Variables.VARIDENT)
print('')
print('Temperature in each altitude level')
print(Variables.XA)
print('')
VARIDENT
[[0 0 0]]
Temperature in each altitude level
[328.16 305.63 284.53 264.74 246.16 228.71 212.31 196.9 182.41 168.83
156.1 144.21 133.41 124.3 117.86 113.89 110.75 108.25 106.26 104.68
103.4 104.5 107.3 110.1 113.5 117. 120.5 125. 128.5 132.4
136. 139.8 143.4 146.2 148.4 150.3 151.8 153. 154.1 155.1
156. 156.83 157.66 158.49 159.32 160.15 160.98 161.81 162.64 163.47
164.3 165.14 165.98 166.82 167.66 168.5 169.34 170.18 171.02 171.86
172.7 173.53 174.36 175.19 176.02 176.85 177.68 178.51 179.34 180.17
181. 185. 189. 193. 197. 201. 208.13 215.27 222.4 229.53
236.67]
Spectroscopy
In this example, we do not include the k-tables in the current working directory because of storage limitations, although these are available upon request.
Here, we include several radiatively active gases: C\(_2\)H\(_2\), C\(_2\)H\(_6\), \(^{12}\)CH\(_4\), \(^{13}\)CH\(_4\), \(^{12}\)CH\(_3\)D, PH\(_3\) and NH\(_3\).
The FWHM of the instrument function was included in the calculation of the k-tables, being 2.5 cm\(^{-1}\).
[6]:
Spectroscopy.summary_info()
INFO :: summary_info :: Spectroscopy_0.py-357 ::
#===== SUMMARY =====#
Spectroscopy_0 instance at memory location 135726490835936
Calculation type ILBL :: (<SpectralCalculationModeEnum.K_TABLES: 0>, ' (k-distribution)')
Number of radiatively-active gaseous species :: 7
Gaseous species :: ['C2H2', 'C2H6', 'CH4 (1)', 'CH4 (2)', 'CH4 (3)', 'PH3', 'NH3']
Number of g-ordinates :: 20
Number of spectral points :: 1496
Wavelength range :: (np.float64(5.0), '-', np.float64(1500.0))
Step size :: 1.0
Spectral resolution of the k-tables (FWHM) :: 2.5
Number of temperature levels :: 15
Temperature range :: (np.float32(70.0), '-', np.float32(300.0))
Number of pressure levels :: 20
Pressure range :: (np.float32(3.0590232e-07), '-', np.float32(7.3890543))
#===================#
2. Running a foward model
After having read the input files and explored some of the most relevant parameters, we can easily perform a forward model.There are two main ways of doing this:
Since all NEMESIS input files are present in the working directory, the forward model could be performed using the nemesisPY.py program from the command window. This is the main program for running a general retrieval, but if NITER = -1 in the .inp file, then a forward model is calculated, not a retrieval. After running the program, the forward model will be stored in the .mre file.
The forward model can be calculated in python using the Forward Model class. In particular, we will use the nemesisfm(), which is the most general version of the forward model. After running this function, the forward model will be stored as an array.
Here, we use the second approach.
[7]:
ForwardModel = ans.ForwardModel_0(runname=runname, Atmosphere=Atmosphere,Surface=Surface,Measurement=Measurement,Spectroscopy=Spectroscopy,Stellar=Stellar,Scatter=Scatter,CIA=CIA,Layer=Layer,Variables=Variables)
SPECONV = ForwardModel.nemesisfm()
INFO :: __init__ :: ForwardModel_0.py-284 :: Checking atmospheric gasses have spectroscopy data.
WARNING :: __init__ :: ForwardModel_0.py-331 :: Not all atmospheric gasses have spectroscopy data.
# WARNING #########################################################################
The following atmospheric gasses ARE NOT PRESENT in the spectroscopy data and WILL NOT CONTRIBUTE TO OPACITY:
C2H4 (id 32) isotopologue 0
C4H2 (id 30) isotopologue 0
H2 (id 39) isotopologue 0
He (id 40) isotopologue 0
To deactivate this warning place a path to a k-table file for these gasses in one of the following locations (depending upon your input file type):
[HDF5 Input]
In the "cirstest.h5" file, add an entry to "/Spectroscopy/LOCATION"
and update "/Spectroscopy/NGAS" appropriately.
[LEGACY Input]
Add an entry to the "cirstest.kls" file.
# END WARNING #####################################################################
INFO :: read_tables :: Spectroscopy_0.py-831 :: Reading table self.LOCATION=PathRedirectList(['../../../archnemesis/Data/reference_tables/ktab_Jupiter_CIRS_nadir/c2h2ZERO_5-1500_2pt5.kta', '../../../archnemesis/Data/reference_tables/ktab_Jupiter_CIRS_nadir/c2h6ZERO_5-1500_2pt5.kta', '../../../archnemesis/Data/reference_tables/ktab_Jupiter_CIRS_nadir/ch4ONE_5-1500_2pt5.kta', '../../../archnemesis/Data/reference_tables/ktab_Jupiter_CIRS_nadir/ch4TWO_5-1500_2pt5.kta', '../../../archnemesis/Data/reference_tables/ktab_Jupiter_CIRS_nadir/ch4THREE_5-1500_2pt5.kta', '../../../archnemesis/Data/reference_tables/ktab_Jupiter_CIRS_nadir/ph3ZERO_5-1500_2pt5.kta', '../../../archnemesis/Data/reference_tables/ktab_Jupiter_CIRS_nadir/nh3ZERO_5-1500_2pt5.kta'], redirects = {}) wavemin=np.float64(5.0) wavemax=np.float64(1399.0)
WARNING :: layer_split :: Layer_0.py-1450 :: from layer_split() :: LAYHT < H(0). Resetting LAYHT
WARNING :: __init__ :: AtmCalc_0.py-254 :: in AtmCalc_0.py file :: THERM requires BINBB disabled - resetting
INFO :: calculate_vertical_cia_opacity :: ForwardModel_0.py-3812 :: Calculating self.CIAX opacity
INFO :: calculate_layer_opacity :: ForwardModel_0.py-3881 :: CIRSrad :: Aerosol optical depths at (np.float64(5.0), ' :: ', array([0.]))
INFO :: calculate_layer_opacity :: ForwardModel_0.py-3899 :: Calculating TOTAL opacity
INFO :: calculate_layer_opacity :: ForwardModel_0.py-3916 :: CIRSradg :: Calculating TOTAL line-of-sight opacity
INFO :: CIRSrad :: ForwardModel_0.py-4379 :: CIRSrad :: IMODM = <PathCalcEnum.THERMAL_EMISSION: 64>
INFO :: calculate_thermal_emission_spectrum :: ForwardModel_0.py-4060 :: CIRSradg :: Calculating THERMAL_EMISSION
[8]:
fig,ax1 = plt.subplots(1,1,figsize=(8,3))
ax1.plot(Measurement.VCONV[:,0],SPECONV[:,0])
ax1.set_xlabel('Wavenumber (cm$^{-1}$)')
ax1.set_ylabel('Radiance (W cm$^{-2}$ sr$^{-1}$ (cm$^{-1}$)$^{-1}$)')
ax1.set_facecolor('lightgray')
ax1.grid()
plt.tight_layout()
Validating the forward model against NEMESIS
In order to test the radiative transfer calculations in archNEMESIS, we compare the computed forward model with the analogous from NEMESIS (i.e., Fortran version).
Because of historical reasons, the values of the spectra in the .mre file are multiplied by 10\(^9\) (i.e., nW cm\(^{-2}\) sr\(^{-1}\) (cm\(^{-1}\))\(^{-1}\)), which needs to be accounted for when comparing the values with the forward model computed with the nemesisfm() function in NemesisPy.
This comparison shows that both versions of the NEMESIS code agree to less than 0.5%.
[9]:
#Reading .mre file computed with the fortran version
lat,lon,ngeom,ny,wave1,specret1,specmeas,specerrmeas,nx,Var,aprprof,aprerr,retprof,reterr = ans.read_mre('cirstest_fortran')
specret1 = specret1 * 1.0e-9
#Making summary plot
fig,(ax1,ax2) = plt.subplots(2,1,figsize=(8,6),sharex=True)
ax1.plot(wave1,specret1,label='NEMESIS',linewidth=2.)
ax1.plot(Measurement.VCONV[:,0],SPECONV[:,0],label='archNEMESIS',c='tab:red',linewidth=1)
ax2.plot(wave1,(specret1[:,0]-SPECONV[:,0])/SPECONV[:,0]*100.,c='black',linewidth=1.)
ax2.set_ylim(-2.,2.)
ax1.grid()
ax2.grid()
ax1.legend()
ax1.set_facecolor('lightgray')
ax2.set_facecolor('lightgray')
ax2.set_xlabel('Wavelength ($\mu$m)')
ax1.set_ylabel('Radiance (W cm$^{-2}$ sr$^{-1}$ (cm$^{-1}$)$^{-1}$)')
ax2.set_ylabel('Difference with Fortran (%)')
plt.tight_layout()
<>:17: SyntaxWarning: invalid escape sequence '\m'
<>:17: SyntaxWarning: invalid escape sequence '\m'
/tmp/ipykernel_2424644/3287963818.py:17: SyntaxWarning: invalid escape sequence '\m'
ax2.set_xlabel('Wavelength ($\mu$m)')
Calculating the jacobian matrix
Similarly, apart from the forward model, archNEMESIS allows the calculation of the jacobian matrix, which includes information about the sensitivity of the spectrum to each of the parameters in the state vector. These calculations are also easily performed using the Forward Model class using the jacobian_nemesis() function. In particular, all the parameters included in the state vector for this particular example can be calculated analytically, although in some cases the derivatives have to be calculated numerically. The output of this function is the modelled measurement vector (i.e., modelled spectrum, since there is only one geometry in this example), and the jacobian matrix (i.e., derivative of the modelled spectrum with respect to the elements of the state vector, which represent the temperature at different altitude/pressure levels).
[10]:
YN,KK = ForwardModel.jacobian_nemesis()
INFO :: jacobian_nemesis :: ForwardModel_0.py-2223 :: Calculating analytical part of the Jacobian :: Calling nemesisfmg
INFO :: read_tables :: Spectroscopy_0.py-831 :: Reading table self.LOCATION=PathRedirectList(['../../../archnemesis/Data/reference_tables/ktab_Jupiter_CIRS_nadir/c2h2ZERO_5-1500_2pt5.kta', '../../../archnemesis/Data/reference_tables/ktab_Jupiter_CIRS_nadir/c2h6ZERO_5-1500_2pt5.kta', '../../../archnemesis/Data/reference_tables/ktab_Jupiter_CIRS_nadir/ch4ONE_5-1500_2pt5.kta', '../../../archnemesis/Data/reference_tables/ktab_Jupiter_CIRS_nadir/ch4TWO_5-1500_2pt5.kta', '../../../archnemesis/Data/reference_tables/ktab_Jupiter_CIRS_nadir/ch4THREE_5-1500_2pt5.kta', '../../../archnemesis/Data/reference_tables/ktab_Jupiter_CIRS_nadir/ph3ZERO_5-1500_2pt5.kta', '../../../archnemesis/Data/reference_tables/ktab_Jupiter_CIRS_nadir/nh3ZERO_5-1500_2pt5.kta'], redirects = {}) wavemin=np.float64(5.0) wavemax=np.float64(1399.0)
WARNING :: layer_split :: Layer_0.py-1450 :: from layer_split() :: LAYHT < H(0). Resetting LAYHT
WARNING :: __init__ :: AtmCalc_0.py-254 :: in AtmCalc_0.py file :: THERM requires BINBB disabled - resetting
INFO :: calculate_vertical_cia_opacity :: ForwardModel_0.py-3812 :: Calculating self.CIAX opacity
INFO :: calculate_layer_opacity :: ForwardModel_0.py-3881 :: CIRSrad :: Aerosol optical depths at (np.float64(5.0), ' :: ', array([0.]))
INFO :: calculate_layer_opacity :: ForwardModel_0.py-3899 :: Calculating TOTAL opacity
INFO :: calculate_layer_opacity :: ForwardModel_0.py-3916 :: CIRSradg :: Calculating TOTAL line-of-sight opacity
INFO :: CIRSrad :: ForwardModel_0.py-4379 :: CIRSrad :: IMODM = <PathCalcEnum.THERMAL_EMISSION: 64>
INFO :: calculate_thermal_emission_spectrum :: ForwardModel_0.py-4060 :: CIRSradg :: Calculating THERMAL_EMISSION
INFO :: nemesisfmg :: ForwardModel_0.py-675 :: Mapping gradients from Layer to Profile
INFO :: nemesisfmg :: ForwardModel_0.py-689 :: Mapping gradients from Profile to State Vector
[11]:
#Plotting the Jacobian matrix
from mpl_toolkits.axes_grid1 import make_axes_locatable
fig,ax1 = plt.subplots(1,1,figsize=(8,3))
im1 = ax1.contourf(Measurement.VCONV[:,0],Atmosphere.P/1.0e5,np.transpose(KK),levels=51,cmap='inferno')
ax1.set_ylim(Atmosphere.P.max()/1.0e5,Atmosphere.P.min()/1.0e5)
ax1.set_xlabel('Wavenumber (cm$^{-1}$)')
ax1.set_ylabel('Pressure (bar)')
ax1.set_yscale('log')
# create an axes on the right side of ax. The width of cax will be 5%
# of ax and the padding between cax and ax will be fixed at 0.05 inch.
divider = make_axes_locatable(ax1)
cax = divider.append_axes("right", size="5%", pad=0.05)
cbar = plt.colorbar(im1, cax=cax)
cbar.set_label('Jacobian matrix (dR/dT)')
[12]:
#Making summary plot
fig,(ax1,ax2) = plt.subplots(2,1,figsize=(8,6),sharex=True)
ax1.plot(wave1,specret1,label='NEMESIS',linewidth=2.)
ax1.plot(Measurement.VCONV[:,0],YN,label='archNEMESIS',c='tab:red',linewidth=1)
ax2.plot(wave1,(specret1[:,0]-SPECONV[:,0])/SPECONV[:,0]*100.,c='black',linewidth=1.)
ax2.set_ylim(-2.,2.)
ax1.grid()
ax2.grid()
ax1.legend()
ax1.set_facecolor('lightgray')
ax2.set_facecolor('lightgray')
ax2.set_xlabel('Wavelength ($\mu$m)')
ax1.set_ylabel('Radiance (W cm$^{-2}$ sr$^{-1}$ (cm$^{-1}$)$^{-1}$)')
ax2.set_ylabel('Difference with Fortran (%)')
plt.tight_layout()
<>:13: SyntaxWarning: invalid escape sequence '\m'
<>:13: SyntaxWarning: invalid escape sequence '\m'
/tmp/ipykernel_2424644/2005142950.py:13: SyntaxWarning: invalid escape sequence '\m'
ax2.set_xlabel('Wavelength ($\mu$m)')
[ ]: