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.
This notebook demonstrates the use of the code used to read FLIR Inc ptw files, calculating the instrument function and applying the instrument function to an example measurement.
In [1]:
from IPython.display import display
from IPython.display import Image
from IPython.display import HTML
import numpy as np
import pandas as pd
import pyradi.ryptw as ryptw
import pyradi.ryplot as ryplot
import pyradi.ryfiles as ryfiles
import pyradi.ryplanck as ryplanck
import matplotlib as mpl
mpl.rc("savefig", dpi=150)
mpl.rc('figure', figsize=(10,8))
%matplotlib inline
pd.set_option('display.width', 150)
In [ ]:
First download the sample file. You have to be connected to the internet to download the file (875 KB). Alternatively copy the file (pyradiSamplePtw.tgz) from the pyradi data directory to the notebook's working directory. After downloading and opening the file there should be three ptw files, an XML file with calibration data, and a few spectral data files in the current working directory.
In [2]:
tgzFilename = 'pyradiSamplePtw.tgz'
destinationDir = '.'
tarFilename = 'pyradiSamplePtw.tar'
url = 'https://raw.githubusercontent.com/NelisW/pyradi/master/pyradi/data/'
dlNames = ryfiles.downloadUntar(tgzFilename, url, destinationDir, tarFilename)
print('filesAvailable are {}'.format(dlNames))
Open one of the ptw files and print the header contents.
In [3]:
# first read the ptw file
ptwfile = './PyradiSampleLWIR.ptw'
outfilename = 'PyradiSampleLWIR.txt'
header = ryptw.readPTWHeader(ptwfile)
ryptw.showHeader(header)
Next load a number of frames from the file: select which frames by listing the required frames (zero-based) in the list of frames:
In [4]:
#loading sequence of frames
framesToLoad = [3,4,10]
data, fheaders = ryptw.getPTWFrames (header, framesToLoad)
print(data.shape)
print(ptwfile)
print('(1) Data shape={}'.format(data.shape))
ifrm = 0
print('Frame {} summary data:'.format(ifrm))
print('Image data:\n{}'.format(data[ifrm]))
print('Sensor Temperature4 {} K '.format(fheaders[ifrm].h_sensorTemp4))
print('FPA Temperature {} K '.format(fheaders[ifrm].h_detectorTemp))
print('Time {}:{}:{} '.format(fheaders[ifrm].h_frameHour,
fheaders[ifrm].h_frameMinute,fheaders[ifrm].h_frameSecond))
The following code attempts to load a frame that does not exist in the input file:
In [5]:
#loading sequence of frames, with an error in request
framesToLoad = [0,4,10]
data, fheaders = ryptw.getPTWFrames (header, framesToLoad)
print(data.shape)
Load single frames and concatenate into a numpy array. The first index into the array is the frame dimension and the second and third indices are row and columns in the image. First load the 0'th element to create the data array and then append the additional frames to the original data array.
In [6]:
# loading single frames
framesToLoad = list(range(1, 101, 1))
frames = len(framesToLoad)
data, fheaders = ryptw.getPTWFrame (header, framesToLoad[0])
for frame in framesToLoad[1:]:
f, fheaders = ryptw.getPTWFrame (header, frame)
data = np.concatenate((data, f))
rows = header.h_Rows
cols = header.h_Cols
img = data.reshape(frames, rows ,cols)
print(img.shape)
In [7]:
display(Image(filename='images/11-radianceEq.PNG'))
This information is summarised from (Sec 9.8) of my book.
The imaging radiometers and spectrometer were calibrated (or characterized) against laboratory blackbody sources. In most cases the calibration is an elaborate process covering several instrument settings, filters, and environmental conditions, but it can be summarized as follows. Calibration entails measuring the instrument output (voltage or digital levels) versus source temperature for a series of source temperatures (ideally over the full dynamic range of the instrument). The set of measurements form a temperature calibration curve (see figure below). When using the instrument, the curve is read 'backward' such that a given instrument voltage returns the appropriate source temperature. If the test object's emissivity is the same as the laboratory source emissivity, the object's temperature will be equal to the source temperature.
Accurate measurement work requires that the spectral sensor response must be known. In this case, the calibration source temperature, together with the spectral sensor response, can be used to calculate the inband source radiance. A new curve is now constructed to relate signal voltage with inband radiance for an extended source completely filling the pixels. This curve can now be used to determine the radiance field incident on the instrument during a measurement. This radiance value lends itself to further analysis.
The figure below shows typical calibration data. The top curve relates source temperature and instrument voltage (digital level). The bottom curve relates source irradiance on the sensor entrance aperture to instrument voltage. For an extended source, the source apparent radiance is related to irradiance by $E=L\omega$, where $\omega$ is the pixel FOV. In this case the instrument was calibrated for three different gain settings (integration times of 30, 120, and 500 $\mu$s) and two different sensor internal temperature conditions (16 and 42.6 $^\circ$C). Observe the importance of calibration at different internal temperatures: the curves (for the same instrument and settings) show significantly different responses. A good strategy to allow for temperature changes in the sensor is to measure the calibration curve at several different internal temperatures and then interpolate between these, according to the actual instrument temperature during the measurement.
The effect of hot optics is also shown in this figure. It is evident that the hot optics flux sets an asymptotically lower measurable flux limit (called the 'floor'). When measuring low-temperature test samples, a small variation in internal temperature shifts the floor up or down, playing havoc with calibration. In this particular instrument setting, an ND2 neutral density filter (0.01 transmittance) was used. Much of the hot optics flux emanates from this filter. So in all fairness, there is little point in using an ND2 filter when measuring a low-temperature target. The situation does arise, however, if there is a requirement to measure both hot and cold test samples without changing filters.
In [9]:
display(Image(filename='images/radiometry08.PNG'))
Load the calibration data, presumably done with a blackbody and at short range.
In [10]:
calData = ryptw.JadeCalibrationData('./lwir100mm150us10ND20090910.xml', '.')
The calibration file contains all the information required to calculate the instrument function of the camera. It lists the filenames for the various spectral parameters (sensor response, optics transmittance, source emissivity, and atmospheric transmittance). The file also contains key sensor parameters such as detector pitch and fill factor and focal length. The most important information is the list of camera digital levels versus source temperature, for one or two instrument temperatures. By convention, the temperature in this file must be in degrees celcius. An example file contents is shown below.
Two elements in this file are not measured but must be determined as part of the data analysis process:
JadeCalibration/Caldatas/Caldata@DlFloor
JadeCalibration/Caldatas/Caldata@Power<JadeCalibration Name="lwir100mm0150us10ND20090910" ID="lwir100mm10ND" Version="" Summary="" >
<SensorResponse Filename="LWIRsensor.txt"/>
<OpticsTransmittance Filename="LW100mmLens.txt"/>
<Filter Filename="LWND10.txt"/>
<SourceEmis Filename="Unity.txt"/>
<AtmoTau Filename="Unity.txt"/>
<DetectorPitch Value="30e-6"/>
<FillFactor Value="0.7"/>
<Focallength Value="0.1"/>
<IntegrationTime Value="150e-6"/>
<Fnumber Value="2"/>
<Nu Min="700" Max="1665" Inc="5" />
<Caldatas>
<Caldata InstrTemperature="17.1" DlFloor="3625" Power="10" >
<CalPoint Temp="50" DL="4571" Comment="temp in deg C" />
<CalPoint Temp="100" DL="5132" />
<CalPoint Temp="150" DL="5906" />
<CalPoint Temp="200" DL="6887" />
<CalPoint Temp="250" DL="8034" />
<CalPoint Temp="300" DL="9338" />
<CalPoint Temp="350" DL="10834" />
<CalPoint Temp="400" DL="12386" />
<CalPoint Temp="450" DL="14042" />
<!-- <CalPoint Temp="500" DL="15324" /> -->
</Caldata>
<Caldata InstrTemperature="34.4" DlFloor="4210" Power="10" >
<CalPoint Temp="50" DL="5477" Comment="temp in deg C" />
<CalPoint Temp="100" DL="6050" />
<CalPoint Temp="150" DL="6817" />
<CalPoint Temp="200" DL="7789" />
<CalPoint Temp="250" DL="8922" />
<CalPoint Temp="300" DL="10262" />
<CalPoint Temp="350" DL="11694" />
<CalPoint Temp="400" DL="13299" />
<CalPoint Temp="450" DL="14921" />
</Caldata>
</Caldatas>
</JadeCalibration>The init_() function loads the data and calculates the mapping functions between digital level, radiance and temperature. Using the spectral curves and DL vs. temperature calibration inputs it calculate the various mapping functions between digital level, radiance and temperature.
When loaded, plot the spectral data for inspection. The various plot functions plot the data to a series of files with hard-coded filenames, built up from the name of the calibration file, with an appropriate added string.
In [11]:
print(calData.Info())
In [12]:
calData.LU.PlotSpectrals()
# HTML('<img src="images/lwir100mm0150us10ND20090910-spectrals.png" width=800 />')
# display(Image(filename='lwir100mm0150us10ND20090910-spectrals.png'))
The following graph plots the source radiance versus source temperature, for the instrument spectral response, as defined in the input files.
In [13]:
calData.LU.PlotCalTempRadiance()
# HTML('<img src="images/lwir100mm0150us10ND20090910-CalTempRadiance.png" width=800 />')
# display(Image(filename='lwir100mm0150us10ND20090910-CalTempRadiance.png'))
The radiance caused by the hot optics in the spectral band is as follows:
In [14]:
calData.LU.PlotCalTintRad()
# HTML('<img src="images/lwir100mm0150us10ND20090910-CalInternal.png" width=800/>')
# display(Image(filename='lwir100mm0150us10ND20090910-CalInternal.png'))
Plot the spectral radiance data for the temperatures used in the calibration process:
In [15]:
calData.LU.PlotCalSpecRadiance()
# HTML('<img src="images/lwir100mm0150us10ND20090910-CalRadiance.png" width=800/>')
# display(Image(filename='lwir100mm0150us10ND20090910-CalRadiance.png'))
The instrument function is nominally a straight line, but flattens out at the low and high ends. At the low end, the (constant) radiance in the optics and instrument dominates the incoming flux on the detector, resulting in an asymptotic lower limit in digital level. At the top end, the detector or digitiser saturates, resulting in an asymptotic line at the maximum recordable level. In practice, the top-end asymptotic line is avoided simply by not using the instrument at these signal levels. The low-end asymptotic line is always present and must always be accounted for.
The low-end asymptotic line is defined by a set of two parameters in the input file (one set each for each instrument temperature):
JadeCalibration/Caldatas/Caldata@DlFloor
JadeCalibration/Caldatas/Caldata@Power
These two parameters are not directly measurable and must be determined by an iterative process. This is done by entering estimates into the data file and then plotting the two graphs below, changing the values until you are satisfied with the graph. The DlFloor parameter determines the asymptotic limit, and the Power value determines the shape of the transition between the asymptotic line and the linear portion of the instrument function.
In [16]:
calData.LU.PlotCalDLRadiance()
# HTML('<img src="images/lwir100mm0150us10ND20090910-CaldlRadiance.png" width=800/>')
# display(Image(filename='lwir100mm0150us10ND20090910-CaldlRadiance.png'))
In [17]:
calData.LU.PlotCalDLTemp()
# HTML('<img src="images/lwir100mm0150us10ND20090910-CalDLTemp.png" width=800 />')
# display(Image(filename='lwir100mm0150us10ND20090910-CalDLTemp.png'))
In [18]:
calData.LU.PlotCalTempDL()
# HTML('<img src="images/lwir100mm0150us10ND20090910-CalTempDL.png" width=800 />')
# display(Image(filename='lwir100mm0150us10ND20090910-CalDLTemp.png'))
Finally, print out the calibration data and instrument function:
In [19]:
print(calData.LU.Info())
The two instrument functions (at low and high instrument temperatures) are known, the instrument function for any arbitrary instrument temperature can now be determined by interpolation (linear in this case).
The following code is used to verify digital-level-values for temperatures at the measured instrument temperatures and then demonstrates the interpolation of the instrument function at an intermediate instrument temperature.
In [20]:
for Tint in [17.1]:
for DL in [4571, 5132, 5906, 6887, 8034, 9338, 10834, 12386, 14042 ]:
print('Tint={} DL={} T={:.1f} L={:.0f}'.format(Tint, DL, calData.LU.LookupDLTemp(DL, Tint)-273.15, calData.LU.LookupDLRad(DL, Tint)))
print('----------------------------------------- ')
for Tint in [34.4]:
for DL in [5477, 6050, 6817, 7789, 8922, 10262, 11694, 13299, 14921 ]:
print('Tint={} DL={} T={:.1f} L={:.0f}'.format(Tint, DL, calData.LU.LookupDLTemp(DL, Tint)-273.15, calData.LU.LookupDLRad(DL, Tint)))
print('----------------------------------------- ')
for Tint in [25]:
for DL in [5477, 6050, 6817, 7789, 8922, 10262, 11694, 13299, 14921 ]:
print('Tint={} DL={} T={:.1f} L={:.0f}'.format(Tint, DL, calData.LU.LookupDLTemp(DL, Tint)-273.15, calData.LU.LookupDLRad(DL, Tint)))
The following measurement observed a large-area blackbody in the lab. The blackbody set point was 150 $^\circ$C and a recording made. The Cedip Altair s/w calculated 157 $^\circ$C in the centre area of the source. The ryptw module calculates a temperature of 149 $^\circ$C.
In this analysis we use the calibration file that was analysed above. The instrument (housing) temperature for the relevant frame is extracted from the header frame. Once the instrument internal temperature is known, it is used to interpolate a curve between the two extreme instrument functions. With the instrument function known at the operational temperature, the source temperature is calculated from the digital levels in the ptw file.
The ptw file may have nill data for the FPA and sensor4 temperatures, from the file description document it appears that a nill value means that it is 'disabled' (what ever that may mean). The example below has nill values.
In [21]:
import os.path
#load the calibration file and calculate the instrument functon
calData = ryptw.JadeCalibrationData('./lwir100mm150us10ND20090910.xml', '.')
#load the ptw file header
ptwfile = 'LWIR-BBref-150C-150us.ptw'
header = ryptw.readPTWHeader(ptwfile)
# showHeader(header)
dfpts = pd.DataFrame()
#loading sequence of frames
framesToLoad = [1,2]
for frame in framesToLoad:
data,fheader = ryptw.getPTWFrame (header, frame)
strTime ='{}:{}:{}'.format(fheader.h_frameHour,
fheader.h_frameMinute,fheader.h_frameSecond)
columns=['Frame','Time','CamT4 C','CamT1 C','FPA K','Row','Col','DL','T deg C','T K','L W/(m2.sr)']
# print('Sensor Temperature4 {} K '.format(fheader.h_sensorTemp4))
# print('FPA Temperature {} K '.format(fheader.h_detectorTemp))
# print('Time {}:{}:{} '.format(fheader.h_frameHour,
# fheader.h_frameMinute,fheader.h_frameSecond))
#get the internal temperature from the header and use here
tempIm = calData.LU.LookupDLTemp(data, header.h_HousingTemperature1-273.15)
# get a radiance image
radIm = calData.LU.LookupDLRad(data, header.h_HousingTemperature1-273.15)
# print('Instrument temperature during measurement was {:.2f} C'.format(header.h_HousingTemperature1-273.15))
# print('Temperature at ({},{})={} C'.format(160,120,tempIm[160,120]-273.15))
# print('Temperature at ({},{})={} C'.format(140,110,tempIm[140,110]-273.15))
for coords in [np.argwhere(np.min(radIm)==radIm)[0],[160,120],[140,110],np.argwhere(np.max(radIm)==radIm)[0]]:
df2 = pd.DataFrame([[frame,strTime,
round(fheader.h_sensorTemp4-273.15,1),
round(header.h_HousingTemperature1-273.15,1),
round(fheader.h_detectorTemp,1),
coords[0],coords[1],int(data[coords[0],coords[1]]),
round(tempIm[coords[0],coords[1]]-273.15,1),
round(tempIm[coords[0],coords[1]],1),
round(radIm[coords[0],coords[1]],2),
]],
columns=columns)
dfpts = dfpts.append(df2,ignore_index=True)
dfpts = dfpts.append(pd.DataFrame([len(columns)*[np.NaN]],columns=columns),ignore_index=True)
I = ryplot.Plotter(4, 1, 1,'', figsize=(6,6))
I.showImage(1, tempIm, ptitle='{}, frame {}, temperature in K'.format(ptwfile[:-4],frame), titlefsize=7, cbarshow=True, cbarfontsize=7)
# I.saveFig('{}-{}.png'.format(os.path.basename(ptwfile)[:-4],frame))
print(dfpts)
In [22]:
# display(Image(filename='LWIR-BBref-150C-150us-1.png'))
HTML('<img src="images/LWIR-BBref-150C-150us-1.png" width=800 />')
Out[22]:
In [23]:
##
pts = dfpts[['T K','DL']].values.T
calData.LU.PlotCalTempDL(plotPoints=pts)
In [24]:
##
pts = dfpts[['DL','L W/(m2.sr)']].values.T
calData.LU.PlotCalDLRadiance(plotPoints=pts)
To follow.
In [25]:
# 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 [ ]: