- 1 Main principles
- 2 Input data for MCYFS
- 3 Pre-processing
- 3.1 Import
- 3.2 Geometric correction
- 3.3 Shortwave radiometric operations (channels 1, 2, and 3A)
- 3.4 Longwave radiometric operations (channels 3B, 4, and 5)
- 3.5 Masking
- 3.6 Compositing
- 4 Quality control
- 5 Post-processing
- 5.1 Preliminary adaptations
- 5.2 Vegetation indicators
- 5.3 Time domain indicators and operations
- 5.4 Time series analysis
- 5.5 Databases with Regional Unmixed Means
- 6 Products overview
The AVHRR or Advanced Very High Resolution Radiometer is the main sensor on-board of the NOAA satellite platforms, which spin around the earth in a near-polar orbit at a height of about 833 km and a frequency of about 14.1 cycles per day. The AVHRR is an imaging sensor which measures the radiance emitted or reflected by the earth-atmosphere system simultaneously in 5 spectral bands, with a spatial resolution of 1.1 km sub-nadir and a swath width of 2922km so that that the NOAA-AVHRR system registers the full earth surface once every day. It was originally designed for meteorological purposes and to complement the equatorial view of the geostationary GOES-platform. The history and applications of the AVHRR-sensor (Advanced Very High Resolution Radiometer) on board of the series of the US NOAA-satellites have been documented by Cracknell (1997). Since 1978, 15 different NOAA-satellites were launched. After the first version (TIROS-N), they were renamed to NOAAn and numbered from n=6 until 19, and so far only one of them failed (NOAA13). Over the years, they were equipped with three different types of AVHRR-sensors (AVHRR/1 to AVHRR/3). The main characteristics of these platforms and of the spectral bands installed on the different AVHRR-types are summarised in the two tables below. Regardless the sensor type, the number of registered and transmitted spectral bands has always been 5, with as classical sequence (for AVHRR/2): VIS, NIR, MIR, TIR1, TIR2. For the older AVHRR/1, which lacked the TIR2-band, a copy of TIR1 was transmitted as a substitute. For the most recent AVHRR/3 which actually has 6 spectral detectors, only the SWIR-band (3A) is transmitted during the daytime orbital trajectory, while at night the classical MIR (3B) is used.
Based on their daytime equatorial crossover times (EQTAN), the NOAA-AVHRR’s are commonly classified into morning, noon, and evening types. The NOAA-AVHRR system also has some important shortcomings, for instance:
- No registration is made of the platform position in space and of its attitude (yaw, pitch, role), which complicates the geometric rectification of the imagery.
- Although all NOAA’s were designed as sun-synchronous systems (with fixed crossover times), in practice orbital drift is inevitable, leading to a gradual delay (advancement) in the crossover times for the afternoon (morning) overpass. For instance, from 1989 to 1995 NOAA11’s equatorial crossover time steadily shifted from 13h40’ to about 17h00’. This orbital drift impacts the observed reflectance values through a gradual change in solar position at the time of the observation. The emitted brightness temperatures also change according to diurnal temperature evolutions.
- The system contains on-board calibration standards for the MIR and TIR channels, but not for the shortwave bands (RED, NIR, SWIR). Hence, the latter must be calibrated by means of more complex and less precise “vicarious” methods: NASA systematically monitors the digital measurements of “stable” targets (water, ice, desert sands,...) and distributes the adjusted calibration factors via dedicated websites. This monitoring is further hampered by the effects of orbital drift.
- On-board storage facilities are limited and only a subsampled version of the global data, with degraded resolution, can be transmitted to earth (GAC, pixels of 5 x 3 km). The 1 km resolution imagery (LAC) is only received when the satellite is in view of a ground station, so much of this precious information is lost. Finally, this “decentralized approach” has given rise to a plethora of different pre-processing algorithms and distribution formats.
However, these deficiencies are partly compensated by that fact that the system already became operational in the early eighties. And beyond the meteorological context, the synoptic AVHRR imagery soon proved very useful as well for the monitoring of the land surfaces at continental to global scales. This gave rise to a number of innovative data sets with global images composited at fixed intervals (dekad to month), such as AVHRR-Pathfinder II (8 km resolution, years 1981-2001) (Ouadrari et al., 2003), GIMMS (8 km, 1981-present) (Tucker et al.,2005) and EROS (1 km, parts of 1992-1995) (Eidenshink and Faundeen, 1994). Already since the eighties, all this new information certainly stimulated the research on global vegetation dynamics and generated new insights on the distribution of ecosystems and their susceptibility to changes caused by climatic and human factors .
Input data for MCYFS
So far, the MCYFS has only exploited the imagery of NOAA’s with a noon orbit, which means that in practice only NOAA7, 9, 11, 14, 16, 18 and 19 are used. All bands are processed for MCYFS. Data from 1989 to 2000 where directly acquired by local antenna, and further completed with data from the receiving station in Dundee. From 2001 till 2009, the data were acquired from FU Berlin station. Data before 1989 were obtained from the receiving station of Dundee. Since 2009, the images are obtained through EUMETSAT Advanced Retransmission Service (EARS) system. The acquisition procedure of EARS is based on seven ground stations (see figure below) that acquire data of three-minute duration each time that the platform passes over the station. The obtained segments are sent to the computation centre of EUMETSAT. Then, because of the overlap between the segments, each of the received data is divided into one-minute segments and a procedure for selecting the best segment takes place: overlaps are removed and a continuous set of “best quality” segments is produced.
The results are thus sent to the EUMETCast’s servers where they are acquired and composed into regional images. Due to satellite failures, there are 16 dekads lacking in the entire period from July 1981 till present.
|The missing dekads are|
Expressed per year, MMDD (MM: month, DD: day) * 1983: 0211 * 1984: 1211 * 1993: 1111, 1121, 1201, 1211 * 1994: 0621 * 2000: 0611, 0621, 1001, 1111, 1121, 1201, 1211, 1221 * 2005: 0701
Each composite contains only data from one sensor. The following scheme is used for the data selection
A general overview of the work flow in the VITO NOAA-AVHRR chain is given in the figure to the right. This chain is executed in a commercial environment (AppWorx).
NOAA-AVHRR is the most problematic sensor in terms of spatial geo-location (no monitoring of aircraft attitude), spectro-temporal behaviour (platform drift, long history of different sensors, no on-board calibration lamps, etc.) and general management (no centralized processing facilities, many different data providers and formats). A dedicated AVHRR pre-processing chain was developed that handles all these different issues. The following steps are discussed in the sections below:
- Geometric correction
- Shortwave radiometric operations (channels 1, 2, and 3A)
- Longwave radiometric operations (channels 3B, 4, and 5)
An overview of all pre-processing products is summarized in products overview.
The NOAA-chain supports 4 input image formats:
- Level1B – AIJ
- Level1B – KLMN
Some deviating formats are also supported: corrupt SHARP1 images from the old SPACE software (Kerdiles, 1997), different packing of the data, etc. Also, the level1B data format differ among providers. Level1B products from Dundee SRS, FU Berlin and CLASS are supported by the current version of the importer. Quality tests are included addressing (1) gaps in the video data, (2) switches from channels 3A to 3B and (3) corrupt PRT data.
|More information on these quality tests|
Gaps in the video data * The time information of each scan line is screened for the occurrence of gaps. * If erroneous time information is detected, the time information is corrected where possible to identify gaps. * Empty lines are added if a gap occurs. * A user specified number of lines before and after the gap can be removed additionally. In the MARS processing, the previous and following line are removed. 3A/3B switch * The software screens the data for the occurrence of a switch between channel 3A/3B. * If only one line differs from the previous or following lines, that line is replaced with an empty line. * If one file contains multiple switches, the file is considered erroneous and is omitted from further processing. * The image is cut off at the switch and only the first part is retained Corrupt PRT data * The quality of the PRT, internal blackbody temperature and cold space temperature values is evaluated by calculating their average and standard deviation for each channel. If the standard deviation is too large, the data from the nearest group of 4 good lines are used. * If the telemetry data is not present on the scanlines, the coefficients provided in the image are used for calibration.
Besides the imagery itself the following auxiliary data are needed in a configuration file:
- The launch and end dates of sensors
- The agreement between NOAA letters and numbers indicating the satellite ID
- The date of switch between channel 3A and 3B for NOAA-16 and -18
Further documentation can be found in the POD- and KLM-Guides (see Links).
Geometric correction is done in different steps. First, the geolocation of the image is derived from information on the position of the satellite using an orbital model. For NOAA-AVHRR this results in an approximate location with an associated error several times larger than the spatial resolution of the imagery due to the lack of information of the attitude of the satellite (up to NOAA-14). Therefore, a second step is inevitable to position the image more accurately. This is done by chip matching the AVHRR imagery to reference images with accurate geolocation.
The orbital model used is the SGP4-model. The output of running the orbital model consists of 6 images containing the latitude (lat), the longitude (lon), the view zenith angle (vza), the view azimuth angle (vaa), the solar zenith angle (sza) and the solar azimuth angle (saa) for the image pixels respectively. The auxiliary data for the model are the Two-line elements (TLE) obtained from Celestrak (see Links) . The orbital model is documented by Walliéz. The interested reader is also referred to the POD- and KLM-Guides (see Links).
A chip matching procedure was developed for the NOAA-AVHRR pre-processing chain. It uses monthly composites of SPOT-VGT images as a reference. The different steps are:
- Find suitable chip locations from a pre-defined list or from a random generator
- Find the corresponding chip on the reference image and the image to correct and measure their spatial distance.
- When sufficient chips are found, a polynomial is estimated for remapping the image to correct.
Detailed information on the procedure is given in the section Chip matching procedure.
Shortwave radiometric operations (channels 1, 2, and 3A)
The radiometric calibration of the optical channels (channels 1, 2, and 3A) converts the digital numbers to radiances and then to reflectance values. In the case of NOAA-AVHRR, there is no onboard calibration device to monitor the responsiveness of the sensors. The calibration coefficients assessment is therefore confined to vicarious calibration using stable targets (desert, bright clouds, ocean, etc.). Many different calibration coefficients sets have been published in literature (e.g. Rao and Chen 1995, 1996, 1999, Vermote and El Saleous 1996, Mitchell et al. 1996, Che and Price 1992, Teillet and Holben 1994, Los 1998, Loeb 1997, Tahnk and Coackley 2001). Swinnen and Veroustraete (2008) compared a number of sets of these calibration coefficients for NOAA-9, 11 and 14 (AIJ) on a stable target in the Namib Desert, and the choice of calibration coefficients for these sensors is based on their results.
|More information on the choice of calibration coefficients for NOAA-AIJ|
NOAA Calibration coefficient set 7 Vermote and El Saleous, 1995 9 Vermote and El Saleous, 1995 11 Mitchell, 1999 14 Vermote and El Saleous, 1995
The calibration method differs between the AVHRR/2 and AVHRR/3 sensors. The latter is calibrated with different coefficients for low and high reflectance values, the so-called dual-slope calibration.
For NOAA-KLMN, calibration is based on the coefficients published at NESDIS (see Links).
Equations for the optical calibration
with d = 1 - 0.01672 * cos(0.9856 * (D-4) * PI/180)
NOAA-KLM (dual slope calibration)
RadianceLow reflectance range:
with Ct calculated from
ReflectanceLow reflectance range:
with Ct calculated from:
For the atmospheric correction SMAC v4 (Rahman and Dedieu, 1994) is used. SMAC requires the following auxiliary data:
- Ozone content is obtained from climatology data of the TOMS data. One global set is available for each month. The monthly data are considered to represent day 15 of the month. Resolution of the dataset is 0.25°.
- Water vapour is derived from climatology data. Resolution of the dataset is 0.25° by 0.25°. Monthly values are available and these represent day 15 of the month.
- Tropospheric aerosol is a simple Gaussian curve from north to south. There is no variation in longitude, nor in time. Resolution of the dataset is 0.25°.
- Stratospheric aerosol for the period January 1991 till December 1999 from Sato et al. (1993) is used. It presents monthly values, which vary only with latitude. This accounts for the high aerosol loads after major volcanic eruptions (e.g. El Chicon and Mt. Pinatubo).
- Pressure is calculated from a DEM by means of the formula:
The SMAC-coefficients for NOAA-7, -17 and -18 were calculated by Olivier Hagolle specifically for the MARS-project. These are now freely available in the SMAC-website (see Links). For NOAA-19, the SMAC-coefficients are not yet acquired. Since the spectral response functions are nearly identical to those of NOAA-18 (see KLM-Guide in Links), this set of SMAC coefficients are used to correct imagery from NOAA-19. The input data is interpolated in space to the resolution of the AVHRR images. The nearest input in time is always used. The climatology datasets of ozone and water vapour are from CTIV and are not documented. More information on the atmospheric correction can be found in the following articles: Rahman, H. and G. Dedieu, 1994; Sato et al. 1993. and Berry et al. 2004.
The Normalized Difference Vegetation Index (NDVI) is calculated once all image segments are radiometrically corrected. It is an index based on RED and NIR reflectance to describe photosynthetic activity. It is also used to select the best observation in the compositing procedure. The calculation is done according to the following equation:
Longwave radiometric operations (channels 3B, 4, and 5)
The TOA-radiances are first converted into brightness temperatures (inversion of Planck’s law) and then combined in an advanced split-window approach (Coll & Caselles, 1997) to retrieve the surface temperature of the land pixels (LST). The emissivities are derived from NDVI.
Equations for the longwave radiometric calibration
(D,A,B)from POD-Guide (see Links).
With (b0,1+b1,b2) from KLM-Guide (see Links).
with A,B from KLM-Guide (see Links).
Land Surface Temperature
Land surface temperatures are estimated from the brightness temperatures in two different TIR-bands (eg. B4 and B5 for AVHRR) with the modified split-window method of Coll & Caselles (1997), which requires the input of the total amount of water vapour (same as for SMAC) and of the emissivities in both bands. The latter can be derived per pixel from the NDVI-values (or strongly correlated measures such as fAPAR or fCover). Sea surface temperatures can be assessed with a simplified version of Coll & Caselles where the emissivities are set at 1. The method is explained in detail in Land Surface temperature (LST).
Reflective part of channel 3B
Channel 3B contains an equal amount of reflected and emitted energy. A separate correction is foreseen to retrieve this reflectance.
The equations to derive the reflectance of channel 3B are described below.
Non-clear vs clear pixels detection
A neural network was trained for each of the NOAA-AVHRR sensors based on the following input layers, Top-of-Canopy RED and NIR reflectance, Brightness temperature from channel 4 and 5 and the month in the year. A normal feed forward network with one hidden layer and 5 nodes was used.
Auxiliary data used in the non-clear pixels detection consists of a land sea mask and text files per NOAA platform containing the weights required by the neural network.
For the AIJ-series the method of Gesell (1989) is used to separate pixels detected by the previous mask into cloud and snow pixels.
For the KLM-series, the third channel is for some periods a SWIR channel, and for other periods a MIR channel. Depending on which part of the electromagnetic spectrum is sensed, another method for snow detection is applied. Also here, the method is only applied on the pixels detected in the non-clear pixels mask.
If channel 3 is MIR, snow detection is performed using the method of Gesell (1999).
If it is a SWIR channel, a method based on the NDSII (Normalised Difference Snow and Ice Index) is used (Hall et al., 1995; Xiangming et al., 2001).
This algoritm is explained here
If the NDSII is larger than 0.70 and the NIR-reflectance exceeds 0.375 then the pixel is flagged as snow.
Auxiliary data are a land sea mask (based on GLC2000) and the non-clear pixels mask.
At the end of the pre-processing, all layers are composited to 10-daily composites using a constrained maximum NDVI approach. The constraints used are:
- Cloudy and snow observations are avoided
- Viewing zenith angle
- Solar zenith angle
A Land/sea mask is needed, which is derived from GLC2000 map.
As quality indicators, a number of additional images are created:
- The number of good observations per compositing period.
- The time grid including segment ID.
- The geolocation quality.
The latter is realized by applying a chip matching procedure using fixed geolocation points for which the distortions are calculated. More information is provided on the page Quality of the compositing.
A detailed description of the compositing method can be found in the section 10-daily composites.
During the compositing, the number of frost days are also computed based on the daily LST images.
After the compositing , a new status map is created. More information on the meaning of the bits can be found below.
The VITO NOAA-chain includes a basic quality control (QC) system. In the QC system quicklooks derived from RGB images using a combination of the different bands, overlaid with a vector with the coastlines and the lakes are produced for all segments as well as for the daily composite images. The quicklooks are stored in a database together with other "quality information" as the RMSE, fill level, etc. of the segments. Using the QC system an operator can quickly screen the pre-processing results and set aside bad segments/composites for further analysis so that these can be recomputed afterwards or removed from the archive. Of course not all types of errors can be detected from these quicklooks but it is an easy way to check for large geometric errors (large shifts) and important radiometrical problems (stripes, detector failures,…).
Examples of errors that are detected by the QC system are shown in the figure below.
All Products and algorithms of the post-processing are explained in detail in the concerned section. A short description of the NOAA-AVHRR post-processing products is provided here.
The information contained in the status mask, derived from the 10-daily compositing, is applied on all composite layers, such that the flag information is easily available for further processing. The procedure is described here.
Smoothing reduces the noise in a pixel’s time series and results in time profiles that are more close to reality instead of containing e.g. undetected clouds. The operation is optional and can only be applied on vegetation indicators, because they have a gradual seasonal evolution. Reflectances cannot be smoothed. The method is described in the section Smoothing.
The fraction of Absorbed Photosynthetically Active Radiation is derived from the 10-daily composites of RED and NIR reflectances using the method described by Weiss et al. (2010). The procedure is described in the section fAPAR. The fAPAR is also used in the estimation of the DMP after smoothing.
Dry Matter Productivity (DMP) expresses the primary productivity (expressed in g DM/ha/day) in a 10-daily time step. The method uses smoothed fAPAR and meteo parameters (from Weather Monitoring), and is based on the Monteith equation explained in detail in the section DMP.
Time domain indicators and operations
The smoothed and non-smoothed vegetation indicators (NDVI, fAPAR, DMP and LST) are further aggregated to monthly composites. A different method is used than for the 10-daily composites and is described here.
Differences to previous year
Two difference operators (absolute and relative) are used to calculate the difference between the current composite and the corresponding one from the previous year. The difference operators are discussed in the section Difference Images.
The historical year is a kind of climatology derived from the entire time series of vegetation indicators (NDVI, fAPAR, LST and DMP) derived from NOAA-AVHRR since mid-1981. The historical year is updated every two years. It contains the average value, the standard deviation, minimum , maximum and deciles per 10% per composite. More information is provided in the section Long Term Statistics.
Difference to historical year
The same difference operators are applied between the current composite and the historical year for the vegetation indicators NDVI, fAPAR, DMP and LST. For NDVI, two additional differences are calculated, i.e. the Vegetation Condition Index (VCI) and the Vegetation Productivity Index (VPI). The methods are described in the section Difference Images.
Time series analysis
The purpose of the similarity analysis is to identify the most similar year compared to the current year. This analysis is performed on the time series evolution per pixel since last October or last March. Only pixels from the classes ‘arable land’, ‘pasture’ and ‘rice’ are considered. The similarity analysis is applied on the smoothed time series of NDVI and fAPAR. The method is described in the section Time series analysis.
The purpose of the cluster classification is to group pixels which have a similar evolution during the current season. It is applied on the smoothed time series of NDVI and fAPAR. The method is described in the section Time series analysis.
Databases with Regional Unmixed Means
Some of the above described products are also delivered in the form of Regional Unmixed Means (RUM) which are then added to the database of the Marsop viewer. Briefly described, RUM-values are derived by averaging (part of) the pixels of a certain administrative region which belong to a certain land cover class. In this way, the derived database can be ingested in a GIS-system. The method to derive the RUM-values is described in the section Databases with Regional Unmixed Means (RUM-values).
The table below provides an overview of all 10-daily and monthly composite products that are generated for MCYFS from NOAA-AVHRR input data.
|Notes on the table|
* Original IMGs labelled with suffix A (SWIR TOC-reflectance) are only available when band 3 is switched to the SWIR-mode (B3A). In the MIR-mode (B3B), only the IMGs with suffix B are delivered (brightness temperature BT3 and the TOA-reflectance in the MIR). * For the actual LST, two QLKs are provided: a normal one (v=t) and a special version (v=tt) with a colour legend adapted to lower temperatures (useful in colder periods/regions). * Meaning of the D-suffixes for the “difference” products: * DIFp: a/r = absolute/relative difference to same period in previous year. * DIFh: 0/1 = absolute/relative difference to same period in historical year (LTA). * DIFh: 3 = relative range in historical year. * DIFh: 4 = historical probability.