This software was developed within the context of the EUMETSAT Satellite Application Facility on Numerical Weather Prediction (NWP SAF), under the Cooperation Agreement dated 1st December 2006, between EUMETSAT and the Met Office, UK, by one or more partners within the NWP SAF. The partners in the NWP SAF are the Met Office, ECMWF, KNMI and Météo-France.
© EUMETSAT. 2006, All Rights Reserved.
Please note that this package requires LAPACK to be installed on your system.
Click here for a pdf version of this file (better for printing).
This software package contains the basic routines for producing reconstructed radiances from spectra obtained from advanced infrared sounders.
Reconstructed radiances are radiances that have been intelligently smoothed such that the atmospheric signal is retained while the instrument noise is suppressed. While the information content of the entire spectrum cannot be increased through the reconstruction process, it allows for the efficient compression of the information from the entire spectrum into a reduced number of channels.
Reconstructed radiances (Antonelli et al., 2004) are formed through the evaluation of the amplitudes, p, of the principal components, L, of the observed spectrum. Here, L is the set of Np leading eigenvectors of the covariance matrix of a representative set of thousands of spectra. p is related to y (the noise normalised radiances with the mean radiance subtracted) through
, are then
calculated from:
to a subset of
NR channels, by replacing the first L above
with LNR, those channels will contain all of the
information present in the Np principal components
provided LNR has ≥ Np positive
singular values. The minimum criterion for this is that
NR ≥ Np and in practice this
criterion is usually sufficient.
The supplied package is split into two parts:
The EOFs, L, used in this process are usually produced from radiance (rather than brightness temperature) spectra normalised by the expected (diagonal) noise. The assumed noise is read into the Create_EOFs from a file and then written to the resulting eigenvector file for use by RR_Filter. If one does not wish to noise normalise, one may simply set the values in the assumed noise file to be all unity.
The radiance units used in the input spectra for Create_EOFs must, of course, be consistent with the units used in the noise file and also with the input spectra used by RR_Filter. If one wants to, say, convert IASI apodised radiances to IASI unapodised radiances or convert brightness temperatures to radiances, these transformations must be made before the radiances are presented to the packages.
Two noise files are provided with the package: one derived before launch (IASI_NOISE_8461_PreLaunch.dat) and one from after the instrument became operational (IASI_NOISE_8461_Nov07.dat). The latter is used in the test scripts.
Create_EOFs is a stand-alone program which creates the eigenvectors used by RR_Filter. The eigenvectors are created from a covariance matrix which describes that variabilty of a training set of input spectra that are read in from a supplied file. These training spectra may be either simulated or observed spectra depending on what is required.
There are two strategies for computing EOFs from a number of spectra. The one used here (computation of a covariance matrix from which the EOFs are calculated) or direct calculation from a matrix of each of the spectra via Singlular Value Decomposition. The former has been chosen here to allow large numbers of spectra (many tens of thousands) to be used in the most efficient manner as the singular value decomposition method is limited by the typically available computer memory once tens of thousands of spectra are beig processed.
Create_EOFs will optionally produce the covariance matrix from input spectra; produce the EOFs from a pre-computed covariance matrix; or both. The data flow for the Create_EOFs program is illustrated below:
2.1.2. Files:
Note that some of these files are described further in
Section 2.2.
2.2 Files
2.2.1. Covariance File:
The covariance file is an unformatted binary file (to allow swifter
input/output and because it is not expected that one might want
to inspect it very often).
It actually contains the number of spectra used, n; the
sum of the radiances for each channel:
The true mean and covariance of the input spectra can then be calculated simply while allowing the covariance to be easily updated with additional spectra:
and
.The files supplied with this package are big-endian, so a suitable compiler flag should be used if the machine being used is little-endian.
2.2.2. Radiances File:
The input radiances file is also a binary, direct-access fortran
file (and the supplied example is big-endian). It is possible that the
user will want to change the form of
the input file to match the particular form of the data available.
The supplied file contains simulated IASI radiances in mW/m2/sr/(cm-1) = erg/s/cm2/sr/(cm-1).
2.2.3. Eigenvectors File:
The eigenvectors file is stored as ASCII (with double-precision
(8-byte) reals). It contains all of the
information required to produce reconstructed radiances from an input
spectrum. The entries in this file are as follows:
2.3. Compiling and Running the code:
The program is contained in a single file except for the LAPACK
routines, so it may be compiled with one line provided
LAPACK and BLAS libraries are available. e.g.,
pgf90 -Mbyteswapio -fast -o Create_EOFs Create_EOFs.f90 -llapack -lblasis the command to compile the code with a Portland Group Fortran 90 compiler run on a little-endian machine.
-Mbyteswapio is the flag used to swap from big- to little- endian. This, (or the equivalent for other machines) is needed if one wishes to use the binary data files (which are big-endian) supplied with this package and one is using a little-endian machine. If one is running on a little-endian machine and using a compiler that does not support byte-swapping, one may request a big-endian dataset from the NWPSAF at http://nwpsaf.eu/feedback.html
The program may then be run simply by executing Create_EOFs.
LAPACK (Anderson et al., 2002) and BLAS (Lawson et al., 1979; Dongarra et al., 1988a, 1988b, 1990; Blackford et al., 2002; Dongarra, 2002) are often installed by default on scientific computing systems but if not available already they may by installed from http://www.netlib.org/lapack and http://www.netlib.org/blas. Both are freely-available software packages which may be incorporated into commercial software packages.
Run time and memory requirements vary according to the number of channels being considered. For large numbers of channels, the dominant step in the eigenvector calculation is the conversion of the covariance matrix to tridiagonal form before deriving the eigenvectors. The CPU time for this step scales as the cube of the number of channels used. The following are guideline timings for 1000 and 8461 channels on an IBM Thinkcentre with a 3.40Ghz Intel Pentium 4 CPU:
| 1000 Channels | 8461 Channels | |
| Covariance calculation: Time per 10000 spectra | 1 minute | 74 minutes |
| Eigenvector calculation (1000 EOFs) | 10seconds | 37 minutes |
| Eigenvector calculation (100 EOFs) | 4 seconds | 26 minutes |
| Memory usage | 10Mb | 548Mb |
If one runs on a IBM Cluster 1600 supercomputer, the timings for the 8461 channels case are reduced to 22, 14 and 10 minutes.
Note that as the purpose of reconstructed radiances is to compress the information available in the whole spectrum into a subset of channels, the input channels should comprise a large fraction of the total spectrum. Reasons for excluding channels might be the poor quality of certain individual channels (e.g., the "popping" channels of AIRS) or excessive oversampling of the spectrum produced from an interferometer.
The three subroutines presented here are used to produced filtered radiances from input unfiltered radiances. RR_Filter does this in one step while Calculate_PCScores and PCScores2Spectrum achieve the same result in two steps through the calculation of the principal component scores (which may be stored more efficiently that either the input or output spectra).
A test program, RR_Filter_Test, that calls these subroutines is supplied. This program may be compiled by modifying and executing the file Make_RR_Filter_Test. On running the resultant executable file, RR_Filter_Test, the test program will process 100 test spectra, reporting the QC index for each and outputting the ASCII file RR_Filter_Test.out which contains the means and standard deviations for the input and output spectra (for the 300 output channels specified in the program).
The subroutine arguments are as follows:
| Spectrum_In(NumChans_In) | REAL Array (IN) | The input spectrum |
| Chans_In(NumChans_In) | INTEGER Array (IN) | The channel numbers for the input spectrum. These channels must include all the channels in the Eigenvector file. |
| NumChans_In | INTEGER (IN) | The number of input channels |
| Num_EOFs | INTEGER (IN) | The number of EOFs to be used |
| Chans_Out(NumChans_Out) | INTEGER Array (IN) | The required channel numbers for the output spectrum |
| NumChans_Out | INTEGER (IN) | The required number of output channels |
| Fixed_Channels | LOGICAL (IN) | Set .TRUE. if and only if the input and output channels are fixed. It is recommended for reasons of efficiency that this is set to .TRUE. if at all possible. |
| Last_Call | LOGICAL (IN) | Set to .TRUE. to DEALLOCATE allocated arrays (the spectrum is not processed if LastCall=.TRUE.) |
| Spectrum_Out(NumChans_Out) | REAL Array (OUT) | The output spectrum |
| QC | REAL (OUT) | Quality control index (see below) |
| ErrorCode | INTEGER (OUT) | Subroutine error code (set to zero if the subroutine completed successfully) |
| PC_Scores(Num_EOFs) | REAL Array (OPTIONAL, OUT) | Amplitudes of Principal Components |
| ErrorMatrix (NumChans_Out, NumChans_Out) | REAL Array (OPTIONAL, OUT) | Estimated error covariance matrix (see below) |
The QC index contains the RMS difference between the input and output spectrum. If the noise on the input spectrum and the assumed noise are similar, the QC values should average to somewhat less than unity. If the input spectra are noise-free (e.g., simulations), the QC values should be small (~0.1).
The estimated error covariance matrix, ErrorMatrix, calculates how the (assumed diagonal) instrument noise is transformed by the reconstruction process. It is calculated by: RRR= LNRLTRLTLNR where R and RRR are the error covariance matrices before and after the reconstruction process and LNR and L are the eigenvectors as defined in the introduction. As L is fixed, RRR will not change as long as the channel choice is constant (so need only be evaluated once).
A namelist file, RR_Filter.NL, may be used to specify the Eigenvectors file being used. The namelist name is RRFilter with one variable: EigenvectorsFile. The default value for EigenvectorsFile is Eigenvectors.out.
A test program, RR_Filter_Test, that calls RR_Filter is supplied. This program may be compiled by modifying and executing the file Make_RR_Filter_Test. On running the resultant executable file, RR_Filter_Test, the test program will process 100 test spectra, reporting the QC index for each and outputting the ASCII file RR_Filter_Test.out which contains the means and standard deviations for the input and output spectra (for the 300 output channels specified in the program).
The following figure shows the standard deviations of the noisy-minus-true spectrum (black); the filtered-minus-true spectrum (blue); and the estimated noise on the filtered spectrum (red). Note that this example uses simulated data and also the test spectra are part of the training set.
The subroutine arguments are as follows:
| Spectrum_In(NumChans_In) | REAL Array (IN) | The input spectrum |
| Chans_In(NumChans_In) | INTEGER Array (IN) | The channel numbers for the input spectrum. These channels must include all the channels in the Eigenvector file. |
| NumChans_In | INTEGER (IN) | The number of input channels |
| Num_EOFs | INTEGER (IN) | The number of EOFs to be used |
| Last_Call | LOGICAL (IN) | Set to .TRUE. to DEALLOCATE allocated arrays (the spectrum is not processed if LastCall=.TRUE.) |
| PC_Scores(Num_EOFs) | REAL Array (OUT) | Amplitudes of Principal Components |
| QC | REAL (OUT) | Quality control index (see below) |
| ErrorCode | INTEGER (OUT) | Subroutine error code (set to zero if the subroutine completed successfully) |
The subroutine arguments are as follows:
| PC_Scores(Num_EOFs) | REAL Array (IN) | Amplitudes of Principal Components |
| Num_EOFs | INTEGER (IN) | The number of EOFs to be used |
| Chans_Out(NumChans_Out) | INTEGER Array (IN) | The required channel numbers for the output spectrum |
| NumChans_Out | INTEGER (IN) | The required number of output channels |
| Fixed_Channels | LOGICAL (IN) | Set .TRUE. if and only if the input and output channels are fixed. It is recommended for reasons of efficiency that this is set to .TRUE. if at all possible. |
| Last_Call | LOGICAL (IN) | Set to .TRUE. to DEALLOCATE allocated arrays (the spectrum is not processed if LastCall=.TRUE.) |
| Spectrum_Out(NumChans_Out) | REAL Array (OUT) | The output spectrum |
| ErrorCode | INTEGER (OUT) | Subroutine error code (set to zero if the subroutine completed successfully) |
| ErrorMatrix (NumChans_Out, NumChans_Out) | REAL Array (OPTIONAL, OUT) | Estimated error covariance matrix (see below) |
Two test scripts are provided for testing. SetUp_Test_Quick.sh uses a reduced number of channels, spectra and eigenvectors and runs in about ten seconds. SetUp_Test.sh is a more comprehensive test and takes 2-3 hours.
This section describes the operation of the SetUp_Test_Quick.sh script including a description of the input and output files.
The script compiles the programs, reads in spectra to produce a covariance matrix, and then computes the leading eigenvectors. The script should be executed from the directory in which it resides.
Two "Make" scripts are executed (they are simply shell scripts rather than proper Make files) to compile Create_EOFs (Make_CreateEOFS) and RR_Filter_Test_Quick (Make_RR_Filter_Test_Quick).
A covariance matrix this then constructed from input spectra. This is done through the Create_EOFs program which is controlled by the Create_EOFs_Test_Quick.NL namelist, reproduced below:
&CreateEOFs RadiancesFile = 'data_quick/RADIANCES_100Subset.dat' CovarianceFile = 'data_quick/Cov_Test.out' ChannelFile = 'Chans2Use_Test.dat' MaxSpec=100 AddNoise=F CreateCovariance=T WriteCovariance=T CalculateEOFs=F /
The Create_EOFs program is now run a second time with the Create_EOFs_Test2_Quick.NL namelist:
&CreateEOFs RadiancesFile = 'data_quick/RADIANCES_100Subsetx.dat' CovarianceFile = 'data_quick/Cov_Test.out' EigenvectorsFile = 'data_quick/Eigenvectors_Test.out' ChannelFile = 'Chans2Use_Test.dat' NumEOFs=100 MaxSpec=100 AddNoise=F AddToExistingCovariance=T CreateCovariance=T WriteCovariance=T CalculateEOFs=T /
The final stage of SetUp_Test_Quick.sh is the generation of filtered radiances through the test program RR_Filter_Test_Quick. The control namelist for this program, RR_Filter_Test_Quick.NL only contains the name of the required eigenvector file, viz:
&RRFilter EigenvectorsFile = 'data_quick/Eigenvectors_Test.out' /
The program calculates the noise-filtered radiances from the first 100 spectra in the noisy radiance file. The noise filtered radiances for a subset of 136 channels are calculated. These parameters being hard-coded in the test program. On generating the filtered radiances, a quality control value is returned to standard out for each spectrum, plus its mean value for all spectra. The statistics of the filtering process are output into the file RR_Filter_Test_Quick.out as follows:
The SetUp_Test.sh is very similar to SetUp_Test_Quick.sh except that all channels and 10000 spectra are used in each of the two Create_EOFs runs; 500 EOFs are calculated in the second of these runs and 300 channels are output from RR_Filter_Test. The output from the latter is RR_Filter_Test.out and may be compared to Test_Runs/RR_Filter_Test.out. The eigenvector file for comparison is at Test_Runs/Eigenvectors_Test.sav which is a link to a file in the data for disk usage reasons.
Memory fault:
A memory fault can arise if the called and calling subroutine
do not have matching argument types. To prevent this occuring
interface blocks have been used throughout. This is particularly
important for RR_Filter where the called
subroutine has an optional argument.
Problems are encountered reading binary files:
If problems are encoutered reading binary files, note
that the files supplied with this package are big-endian, so a
suitable compiler flag should be used if the machine being used is
little-endian.
Also, it is assumed that when opening a direct access file (such as the radiances file) the record length "RECL" is given in units of bytes. Some compilers use 4-byte units by default, so you will either have to modify the code or use an appropriate compiler directive (e.g. for ifort you can use the directive "-assume byterecl").
Test runs are not bit identical to supplied results:
Results will not be bit identical to the output in
Test_Runs/ as floating point calculations
depend in detail on the machine, compiler and compiler options used.
Agreement in the reconstructed radiances should be at least 5
sig.figs. and there should be similar accuracy in the leading
eigenvectors.
Package is being run on a little-endian machine and compiler lacks a byte-swapping option: If one is running on a little-endian machine and using a compiler that does not support byte-swapping, one may request a big-endian dataset from the NWPSAF at http://nwpsaf.eu/feedback.html
You can contact the NWPSAF team using the form at: http://nwpsaf.eu/feedback.html
Anderson, E. and Bai, Z. and Bischof, C. and Blackford, S. and Demmel, J. and Dongarra, J. and Du Croz, J. and Greenbaum, A. and Hammarling, S. and McKenney, A. and Sorensen, D., LAPACK Users' Guide, Third Edition, Society for Industrial and Applied Mathematics, Philadelphia, PA, 1999, ISBN: 0-89871-447-8 (paperback)
P. Antonelli, H.E. Revercomb, L.A. Sromovsky, W.L. Smith, R.O. Knuteson, D.C. Tobin, R.K. Garcia, H.B. Howell, H.-L. Huang, and F.A. Best (2004). A principal component noise filter for high spectral resolution infrared measurements, J. Geophys. Res., 109, D23102-23124.
L. S. Blackford, J. Demmel, J. Dongarra, I. Duff, S. Hammarling, G. Henry, M. Heroux, L. Kaufman, A. Lumsdaine, A. Petitet, R. Pozo, K. Remington, R. C. Whaley, An Updated Set of Basic Linear Algebra Subprograms (BLAS), ACM Trans. Math. Soft., 28-2 (2002), pp. 135-151.
J. J. Dongarra, J. Du Croz, S. Hammarling, and R. J. Hanson, An extended set of FORTRAN Basic Linear Algebra Subprograms, ACM Trans. Math. Soft., 14 (1988a), pp. 1-17.
J. J. Dongarra, J. Du Croz, S. Hammarling, and R. J. Hanson, Algorithm 656: An extended set of FORTRAN Basic Linear Algebra Subprograms, ACM Trans. Math. Soft., 14 (1988b), pp. 18-32.
J. J. Dongarra, J. Du Croz, I. S. Duff, and S. Hammarling, A set of Level 3 Basic Linear Algebra Subprograms, ACM Trans. Math. Soft., 16 (1990), pp. 1-17.
J. J. Dongarra, J. Du Croz, I. S. Duff, and S. Hammarling, Algorithm 679: A set of Level 3 Basic Linear Algebra Subprograms, ACM Trans. Math. Soft., 16 (1990), pp. 18-28.
J. Dongarra, Basic Linear Algebra Subprograms Technical Forum Standard, International Journal of High Performance Applications and Supercomputing, 16(1) (2002), pp. 1-111, and International Journal of High Performance Applications and Supercomputing, 16(2) (2002), pp. 115-199.
C. L. Lawson, R. J. Hanson, D. Kincaid, and F. T. Krogh, Basic Linear Algebra Subprograms for FORTRAN usage, ACM Trans. Math. Soft., 5 (1979), pp. 308-323.