This notebook forms part of a series on computational optical radiometry. The notebooks can be downloaded from Github. These notebooks are constantly revised and updated, please revisit from time to time.
The date of this document and module versions used in this document are given at the end of the file.
Feedback is appreciated: neliswillers at gmail dot com.
The pyradi library has a module to handle Modtran files. The module currently has only one function: to read Modtran tape7 files in an intelligent manner. The function is rymodtran.loadtape7.
In [2]:
get_ipython().run_line_magic('matplotlib', 'inline')
import numpy as np
import pandas as pd
import pyradi.ryfiles as ryfiles
import pyradi.ryplot as ryplot
import pyradi.rymodtran as rymodtran
import pyradi.ryutils as ryutils
from IPython.display import display
from IPython.display import Image
from IPython.display import HTML
import matplotlib as mpl
mpl.rc("savefig", dpi=300)
mpl.rc('figure', figsize=(10,8))
# %config InlineBackend.figure_format = 'svg'
# import pandas as pd
# pd.set_option('display.max_columns', 80)
# pd.set_option('display.width', 100)
# pd.set_option('display.max_colwidth', 150)
This function reads in the tape7 file from MODerate spectral resolution atmospheric TRANsmission (MODTRAN) code, that is used to model the propagation of the electromagnetic radiation through the atmosphere. tape7 is a primary file that contains all the spectral results of the MODTRAN run. The header information in the tape7 file contains portions of the tape5 information that is ignored here. The header section in tape7 is followed by a list of spectral samples with corresponding transmissions. Each column has a different component of the transmission or radiance. For more detail, see the Modtran documentation.
The user selects the appropriate columns by listing the column names, as listed below. The tape7.scn file has missing columns, so this function does not work for tape7.scn files. If you need a tape7.scn file with all the columns populated you would have to use the regular tape7 file and convolve this to lower resolution.
The format of the tape7 file changes for different IEMSCT values. For the most part the differences are hidden in the details, if the columns are read by column header and not by counting columns. The various column headers used in the tape7 file are as follows:
IEMSCT = 0 has two lines of column headers. In order to select the column, you must concatenate the two column headers with an underscore in between. All columns are available with the following column names:
FREQ_CM-1, COMBIN_TRANS, H2O_TRANS, UMIX_TRANS, O3_TRANS, TRACE_TRANS,
N2_CONT, H2O_CONT, MOLEC_SCAT, AER+CLD_TRANS, HNO3_TRANS,
AER+CLD_abTRNS, -LOG_COMBIN, CO2_TRANS, CO_TRANS, CH4_TRANS,
N2O_TRANS, O2_TRANS, NH3_TRANS, NO_TRANS, NO2_TRANS,
SO2_TRANS, CLOUD_TRANS, CFC11_TRANS, CFC12_TRANS, CFC13_TRANS,
CFC14_TRANS, CFC22_TRANS, CFC113_TRANS, CFC114_TRANS,
CFC115_TRANS, CLONO2_TRANS, HNO4_TRANS, CHCL2F_TRANS,
CCL4_TRANS, N2O5_TRANS
IEMSCT = 1 has single line of column headers. A number of columns have
headers, but with no column numeric data. In the following list the
columns with header names ** are empty in some runs:
FREQ,TOT_TRANS, PTH_THRML, THRML_SCT, SURF_EMIS, *SOL_SCAT*,
*SING_SCAT*, GRND_RFLT, *DRCT_RFLT*, TOTAL_RAD, *REF_SOL*, *SOL@OBS*,
DEPTH, DIR_EM, *TOA_SUN*, BBODY_T[K]
Hence, these columns may not have valid data:
SOL_SCAT, SING_SCAT, DRCT_RFLT, REF_SOL,SOL@OBS, TOA_SUN
IEMSCT = 2 has a single line column headers. All the columns are available:
FREQ, TOT_TRANS, PTH_THRML, THRML_SCT, SURF_EMIS, SOL_SCAT,
SING_SCAT, GRND_RFLT, DRCT_RFLT, TOTAL_RAD, REF_SOL, SOL@OBS,
DEPTH, DIR_EM, TOA_SUN, BBODY_T[K]
IEMSCT = 3 has a single line column headers. One of these seems to be two words, which, in this code must be concatenated with an underscore. There is also additional column (assumed to be depth in this code). The columns available are
FREQ, TRANS, SOL_TR, SOLAR, DEPTH
First prepare by importing and downloading the tape7 files from the internet.
In [3]:
get_ipython().system('dir')
tgzFilename = 'modtrandata.tgz'
destinationDir = '.'
tarFilename = 'modtrandata.tar'
url = 'https://raw.githubusercontent.com/NelisW/pyradi/master/pyradi/data/'
dlNames = ryfiles.downloadUntar(tgzFilename, url, destinationDir, tarFilename)
print('filesAvailable are {}'.format(dlNames))
The function rymodtran.loadtape7 is relatively easy to use, with only two parameters: the filename and the specification of which columns to read (a list). In some instances (IEMSCT == 0) the column header is given over two lines, and your specification must include the information on both lines. This is done by concatenating the two parts with an underscore character as in TRACE and TRANS becoming TRACE_TRANS.
IEMSCT = 0: this example loads data from five columns. Note the underscore to indicate that the column header flows over two lines.
In [4]:
tape7 = rymodtran.loadtape7("tape7-01", ['FREQ_CM-1', 'COMBIN_TRANS',
'MOLEC_SCAT', 'AER+CLD_TRANS',
'AER+CLD_abTRNS'] )
print(tape7.shape)
The next example reads all the columns in the tape7 file:
In [5]:
tape7 = rymodtran.loadtape7("tape7-01", ['FREQ_CM-1', 'COMBIN_TRANS',
'H2O_TRANS', 'UMIX_TRANS',
'O3_TRANS', 'TRACE_TRANS',
'N2_CONT', 'H2O_CONT', 'MOLEC_SCAT',
'AER+CLD_TRANS', 'HNO3_TRANS',
'AER+CLD_abTRNS', '-LOG_COMBIN',
'CO2_TRANS', 'CO_TRANS', 'CH4_TRANS',
'N2O_TRANS', 'O2_TRANS', 'NH3_TRANS',
'NO_TRANS', 'NO2_TRANS', 'SO2_TRANS',
'CLOUD_TRANS', 'CFC11_TRANS',
'CFC12_TRANS', 'CFC13_TRANS',
'CFC14_TRANS', 'CFC22_TRANS',
'CFC113_TRANS', 'CFC114_TRANS',
'CFC115_TRANS', 'CLONO2_TRANS',
'HNO4_TRANS', 'CHCL2F_TRANS',
'CCL4_TRANS', 'N2O5_TRANS'] )
print(tape7.shape)
IEMSCT = 1: this example loads data from all columns with data in the sample file.
In [6]:
tape7 = rymodtran.loadtape7("tape7-02", ['FREQ', 'TOT_TRANS', 'PTH_THRML',
'THRML_SCT', 'SURF_EMIS',
'GRND_RFLT', 'TOTAL_RAD', 'DEPTH',
'DIR_EM', 'BBODY_T[K]'] )
print(tape7.shape)
IEMSCT = 2: this example loads data from all columns in the sample file.
In [7]:
tape7 = rymodtran.loadtape7("tape7-03", ['FREQ', 'TOT_TRANS', 'PTH_THRML',
'THRML_SCT', 'SURF_EMIS', 'SOL_SCAT',
'SING_SCAT', 'GRND_RFLT',
'DRCT_RFLT', 'TOTAL_RAD', 'REF_SOL',
'SOL@OBS', 'DEPTH', 'DIR_EM',
'TOA_SUN', 'BBODY_T[K]'] )
print(tape7.shape)
IEMSCT = 3: this example loads data from all columns with data in the sample file. Note that in the data file the the fifth column is not labelled, but the function reads still it.
In [8]:
tape7 = rymodtran.loadtape7("tape7-04", ['FREQ', 'TRANS', 'SOL_TR',
'SOLAR', 'DEPTH'] )
print(tape7.shape)
The next example reads data from the file and display the data in a plot. The columns to be read and plotted are defined in the list colSelect. This list is used when loading the data, but also as legend in the plot. The Modtran data is given in wavenumber domain. The plot must be in wavelength domain, hence wavenumber is converted to wavelength using ryutils.convertSpectralDomain. Finally, the various transmittance curves are plotted against wavelength.
In [9]:
colSelect = ['FREQ_CM-1', 'COMBIN_TRANS', 'H2O_TRANS', 'UMIX_TRANS',
'O3_TRANS', 'H2O_CONT', 'MOLEC_SCAT', 'AER+CLD_TRANS']
tape7= rymodtran.loadtape7("tape7VISNIR5kmTrop23Vis", colSelect )
wavelen = ryutils.convertSpectralDomain(tape7[:,0], type='nl')
mT = ryplot.Plotter(1, 1, 1,"Modtran Tropical, 23 km Visibility (Rural)"
+ ", 5 km Path Length",figsize=(12,6))
mT.plot(1, wavelen, tape7[:,1:], "","Wavelength [$\mu$m]", "Transmittance",
label=colSelect[1:],legendAlpha=0.5, pltaxis=[0.4,1, 0, 1]);
The next example plots the transmittance from the individual molecules and aerosol separately to clearly show the effect of each.
In [10]:
# this example plots the individual transmittance components
colSelect = ['FREQ_CM-1', 'COMBIN_TRANS', 'MOLEC_SCAT', 'CO2_TRANS',
'H2O_TRANS', 'H2O_CONT', 'CH4_TRANS', 'O3_TRANS', 'O2_TRANS',
'N2O_TRANS', 'AER+CLD_TRANS', 'SO2_TRANS']
tape7= rymodtran.loadtape7("horizon5kmtropical.fl7", colSelect )
wavelen = ryutils.convertSpectralDomain(tape7[:,0], type='nl')
mT = ryplot.Plotter(1, 9, 1,"Modtran Tropical, 23 km Visibility (Rural)" + ", 5 km Path Length",figsize=(12,21))
mT.semilogX(1, wavelen, tape7[:,1], '','', '', maxNY=4,
label=colSelect[1:],legendAlpha=0.5, pltaxis=[0.2,15, 0, 1]);
mT.semilogX(2, wavelen, tape7[:,2], '','', '', maxNY=4,
label=colSelect[2:],legendAlpha=0.5, pltaxis=[0.2,15, 0, 1]);
mT.semilogX(3, wavelen, tape7[:,10], '','', '', maxNY=4,
label=colSelect[10:],legendAlpha=0.5, pltaxis=[0.2,15, 0, 1]);
mT.semilogX(4, wavelen, tape7[:,4] , '','', '', maxNY=4,
label=colSelect[4:],legendAlpha=0.5, pltaxis=[0.2,15, 0, 1]);
mT.semilogX(5, wavelen, tape7[:,5] , '','', '', maxNY=4,
label=colSelect[5:],legendAlpha=0.5, pltaxis=[0.2,15, 0, 1]);
mT.semilogX(6, wavelen, tape7[:,3] , '','', '', maxNY=4,
label=colSelect[3:],legendAlpha=0.5, pltaxis=[0.2,15, 0, 1]);
mT.semilogX(7, wavelen, tape7[:,6] , '','', '', maxNY=4,
label=colSelect[6:],legendAlpha=0.5, pltaxis=[0.2,15, 0, 1]);
mT.semilogX(8, wavelen, tape7[:,7] * tape7[:,8] , '','', '', maxNY=4,
label=colSelect[7:],legendAlpha=0.5, pltaxis=[0.2,15, 0, 1]);
mT.semilogX(9, wavelen, tape7[:,9] , '','', '', maxNY=4,
label=colSelect[9:],legendAlpha=0.5, pltaxis=[0.2,15, 0, 1]);
In [98]:
# this example plots the individual transmittance components
import matplotlib.patches as patches
colSelect = ['FREQ_CM-1', 'COMBIN_TRANS']
tape7= rymodtran.loadtape7("horizon5kmtropical.fl7", colSelect )
wavelen = ryutils.convertSpectralDomain(tape7[:,0], type='nl')
mT = ryplot.Plotter(1, 1, 1,"Tropical: 27 $^\circ$C, 75% RH, 23 km Visibility (Rural)"
+ ", 5 km Path Length",
figsize=(12,5),doWarning=False)
mT.plot(1, wavelen, tape7[:,1], '','Wavelength [$\mu$m]', 'Transmittance', maxNY=4,maxNX=15,
legendAlpha=0.5, pltaxis=[0.2,15, 0, 1],plotCol=['k']);
ax = mT.getSubPlot(1)
for txt,tp in zip(['NIR','SWIR','MWIR','LWIR'],[(0.04,.048),(0.09,.095),(0.187,.2025),(0.392,.5)]):
p = patches.Rectangle((tp[0],.85),tp[1],.14,fill=True,color='white',transform=ax.transAxes, clip_on=False)
p.set_edgecolor('black')
ax.add_patch(p)
ax.annotate(txt, xy=(tp[0]+tp[1]/2,.9), xycoords='axes fraction', fontsize=16,
horizontalalignment='center', verticalalignment='center')
mT.saveFig('atmosphere2.pdf')
The final example reads path radiance data from a file, plots the results, and calculates the total path radiance in a spectral band. In this example the plot and integration are done in wavenumber domain. Normally you would multiply with a sensor spectral response before integration, but this calculation is over the whole band, equally weighted.
Modtran's path radiance is given in units of W/(cm$^2$.sr) and must ideally be converted to units of W/(m$^2$.sr). This is not a critical conversion, but my preference is to work in base units such as metres.
In [12]:
colSelect = ['FREQ', 'PTH_THRML','SOL_SCAT','SING_SCAT', 'TOTAL_RAD']
skyrad= rymodtran.loadtape7("NIRscat.fl7", colSelect )
sr = ryplot.Plotter(1, 4,1,"Path Radiance in NIR, Path to Space from 3 km",
figsize=(12,8))
# plot the components separately
for i in [1,2,3,4]:
Lpath = 1.0e4 * skyrad[:,i]
sr.plot(i, skyrad[:,0], Lpath, "","Wavenumber [cm$^{-1}$]",
"L [W/(m$^2$.sr.cm$^{-1}$)]",
label=[colSelect[i][:]],legendAlpha=0.5, #pltaxis=[0.4,1, 0, 1],
maxNX=10, maxNY=4, powerLimits = [-4, 4, -5, 5]);
#convert from /cm^2 to /m2 and integrate using the wavenumber vector
totinband = np.trapz(Lpath.reshape(-1, 1),skyrad[:,0],
axis=0)[0]
print('{0} integral is {1} [W/(m^2.sr)]'.format(colSelect[i][:],totinband))
print('Note that multiple scatter contributes significantly' + ' to total path radiance')
Repeat the previous example, but this time plot and integrate in the wavelength domain. The integrated results should be the same as for the example above. In this example the wavenumber domain is converted to wavelength domain and the spectral densities are likewise converted.
Note that the domain and densities are simply converted, not resampled. This means that the data is still at constant wavenumber intervals (varying wavelength intervals). This does not pose a problem when plotting, but care must be taken when integrating to tell trapz that the x interval is not constant. Note, however that the data is in a sequence of increasing wavenumber; that means decreasing wavelength, and as a result the integral is negative (negative integration intervals).
In [14]:
colSelect = ['FREQ', 'PTH_THRML','SOL_SCAT','SING_SCAT', 'TOTAL_RAD']
skyrad= rymodtran.loadtape7("NIRscat.fl7", colSelect )
sr = ryplot.Plotter(1, 4,1,"Path Radiance in NIR, Path to Space from 3 km",
figsize=(12,8))
# plot the components separately
for i in [1,2,3,4]:
wl, Lpath = ryutils.convertSpectralDensity(skyrad[:,0],
skyrad[:,i],'nl')
Lpath *= 1.0e4
sr.plot(i, wl, Lpath, "","Wavelength [$\mu$m]",
"L [W/(m$^2$.sr.$\mu$m)]",
label=[colSelect[i][:]],legendAlpha=0.5, #pltaxis=[0.4,1, 0, 1],
maxNX=10, maxNY=4, powerLimits = [-4, 4, -5, 5],xAxisFmt='%.2f');
#convert from /cm^2 to /m2 and integrate using the wavenumber vector
totinband = - np.trapz(Lpath.reshape(-1, 1),wl,
axis=0)[0]
print('{0} integral is {1} [W/(m^2.sr)]'.format(colSelect[i][:],totinband))
print('Note that multiple scatter contributes significantly' + ' to total path radiance')
Read a file with solar irradiance
In [26]:
colSelect = ['FREQ', 'TRANS','SOL_TR','SOLAR']
sunrad= rymodtran.loadtape7("tape7-sunE", colSelect )
sr = ryplot.Plotter(1,2,1,"Solar Irradiance in NIR", figsize=(12,8))
wl, LsunE = ryutils.convertSpectralDensity(sunrad[:,0],sunrad[:,2],'nl')
wl, LsunS = ryutils.convertSpectralDensity(sunrad[:,0],sunrad[:,3],'nl')
LsunE *= 1.0e4
LsunS *= 1.0e4
sr.plot(1, wl, LsunE, "","Wavelength [$\mu$m]",
"L [W/(m$^2$.sr.$\mu$m)]",
label=['Irradiance on earth'],legendAlpha=0.5, #pltaxis=[0.4,1, 0, 1],
maxNX=10, maxNY=4, powerLimits = [-4, 4, -5, 5],xAxisFmt='%.2f');
sr.plot(1, wl, LsunS, "","Wavelength [$\mu$m]",
"L [W/(m$^2$.sr.$\mu$m)]",
label=['Irradiance outside atmosphere'],legendAlpha=0.5, #pltaxis=[0.4,1, 0, 1],
maxNX=10, maxNY=4, powerLimits = [-4, 4, -5, 5],xAxisFmt='%.2f');
sr.plot(2, wl, sunrad[:,1], "","Wavelength [$\mu$m]",
"$\\tau$",
label=['Transmittance'],legendAlpha=0.5, #pltaxis=[0.4,1, 0, 1],
maxNX=10, maxNY=4, powerLimits = [-4, 4, -5, 5],xAxisFmt='%.2f');
If the spectral transmittance at one range is known, it can be scaled to other distances by applying Beer/Bouguer's law. Calculate the attenuation coefficient for the known transmittance and then apply it to other distances. $$ \tau_a(R_1) = \exp(-R_1 \alpha_a) $$ $$\alpha_a = -\log(\tau_a(R_1))/R_1 $$ and then calculate the transmittance at the required range by $$ \tau_a(R) = \exp(-R \alpha_a) $$
Path radiance can be scaled by multiplication/division by $(1-\tau)$. Divide the known path radiance by $(1-\tau_\textrm{known})$ and then multiply with $(1-\tau_\textrm{new})$. The path radiance for an infinitely long path can be calculated and then rescaled:
$$ L_\textrm{path}(R\rightarrow \infty) = \frac{L_\textrm{path1}}{1-\tau_\textrm{path1}} = \frac{L_\textrm{path2}}{1-\tau_\textrm{path2}} $$Modtran was used to calculate path radiance and transmittance for several distances, using the same atmospheric model. The path is horizontal, at sea level. The Tropical climate model is used with a Rural aerosol and 23 km visibility.
In [24]:
pranges = [100,1000,10000]
p = ryplot.Plotter(1,2,2,figsize=(16,8))
for prange in pranges:
colsPath = ['FREQ', 'TOT_TRANS', 'TOTAL_RAD','DEPTH']
tape7 = rymodtran.loadtape7(f"data/path{prange}/tape7", colsPath)
dfPath = pd.DataFrame(tape7,columns=colsPath)
dfPath = dfPath.rename(columns={'FREQ': 'nu','TOT_TRANS':'tauPath','TOTAL_RAD':'radPathNu'})
wlp, dfPath['radPathWl'] = ryutils.convertSpectralDensity(dfPath['nu'].values,dfPath['radPathNu'].values,'nl')
dfPath['radPathWl'] *= 1.0e4
dfPath['tauPath'] += 1e-6 # prevent zero transmittance
dfPath['wl'] = 1e4 / dfPath['nu']
# dfPath['alphaP'] = -np.log(dfPath['tauPath']) / prange
dfPath['alphaP'] = dfPath['DEPTH'] / prange
dfPath['radPathWlNorm'] = dfPath['radPathWl'] / (1 - dfPath['tauPath'] )
# print(dfPath[dfPath['tauPath']<0.01])
p.plot(1,dfPath['wl'],dfPath['radPathWl'],'Path radiance over distance','Wavelength [$\mu$m]','Radiance W/(m2.sr.um)',
xAxisFmt='%.2f',label=[f'{prange} m'])
p.plot(2,dfPath['wl'],dfPath['tauPath'],'Path transmittance over distance','Wavelength [$\mu$m]','Transmittance',
xAxisFmt='%.2f',label=[f'{prange} m'])
p.plot(3,dfPath['wl'],dfPath['radPathWlNorm'],'Path radiance / (1-$\\tau$)','Wavelength [$\mu$m]','Radiance W/(m2.sr.um)',
xAxisFmt='%.2f',label=[f'{prange} m'])
p.semilogY(4,dfPath['wl'],dfPath['alphaP'],'Path attentuation coefficient','Wavelength [$\mu$m]','Attenuation [1/m]',
xAxisFmt='%.2f',label=[f'{prange} m'])
The bottom-left graph indicates that path radiance can be accurately re-calculated for different distances,
It is evident from the bottom-right graph that the attenuation coefficient approximation is not quite valid to scale for range in strong attenuation lines. The Modtran documentation also mentions somewhere that transmittance at closer ranges is less accurate than at longer ranges.
In spectroscopy the deviations (called nonlinearity errors) from Beer's law is well studied and documented. There are several causes identified most of which is relevant to spectroscopy only (instrument functions, chemical interactions, and what we call path radiance). One cause that could apply here in a modified form is that attenuation decreases at higher concentractions of molecules. At higher concentration (number of molecules along the path?) there is interaction between absorbing particles such that the absorption decreases.
The decreased absorption for longer paths is evident in the bottom-right graph. Note that the difference only occurs in the molecular lines (oxygen in this case). A wavelengths with no molecular absorption, Beer's law holds well. Extending the interaction argument this also means that the effect will be more pronounced for stronger absorbing molecules, such as the oxygen shown here or CO_\textsubscript{2} at 4.3 $\mu$m.
Suppose that the 1000~m transmittance is used as the baseline, then transmittance estimates at shorter distances will be higher than if calculated directly. Conversely transmittance at longer wavelengths will be lower than if calculated directly.
See these and several other sources on the Internet:
https://wikivisually.com/wiki/Beer%E2%80%93Lambert_law
http://simulab.ltt.com.au/5/Laboratory/StudyNotes/snDeviatFromBeerLaw.htm
In [ ]:
In [15]:
# to get software versions
# https://github.com/rasbt/watermark
# you only need to do this once
# pip install watermark
%load_ext watermark
%watermark -v -m -p numpy,scipy,pyradi -g
In [ ]: