Automated Improvement of Geolocation Accuracy in AVHRR Data Using a Two-Step Chip Matching Approach — A Part of the TIMELINE Preprocessor

The geolocation of Advanced Very High Resolution Radiometer (AVHRR) data is known to be imprecise due to minor satellite position and orbit uncertainties. These uncertainties lead to distortions once the data are projected based on the provided orbit parameters. This can cause geolocation errors of up to 10 km per pixel which is an obstacle for applications such as time series analysis, compositing/mosaicking of images, or the combination with other satellite data. Therefore, a fusion of two techniques to match the data in orbit projection has been developed to overcome this limitation, even without the precise knowledge of the orbit parameters. Both techniques attempt to find the best match between small image chips taken from a reference water mask in the first, and from a median Normalized Difference Vegetation Index (NDVI) mask in the second round. This match is determined shifting around the small image chips until the best correlation between reference and satellite data source is found for each respective image part. Only if both attempts result in the same shift in any direction, the position in the orbit is included in a third order polynomial warping process that will ultimately improve the geolocation accuracy of the AVHRR data. The warping updates the latitude and longitude layers and the contents of the data layers itself remain untouched. As such, original sensor measurements are preserved. An included automated quality assessment generates a quality layer that informs about the reliability of the matching.


Introduction
Since 1978, the National Oceanic and Atmospheric Administration (NOAA) operated the Advanced Very High Resolution Radiometer (AVHRR) onboard its Television Infrared Observation Satellite (TIROS) and NOAA satellite platforms.Since 2006, the AVHRR sensors are also mounted aboard the Meteorological Operational Satellites (MetOp) operated by Eumetsat, delivering daily observations of the globe in visible, near infrared, and thermal wavelengths.Three different generations of AVHRR sensors have been developed over the decades, providing observations from four, five, and six spectral channels, respectively.The resolution of the High-Resolution Picture Transmission (HRPT) format data at nadir is around 1.1 km but widens up to 6.5 km with the observation angle, which spans to a maximum of ±55.4 degrees from the nadir [1].Table 1 contains an overview of the different platforms, sensor characteristics, and typical use cases for the available bands.AVHRR data come with calibration coefficients and orbital-elements ("two-line-elements", TLE), which base on the satellite ephemeris model and the instrument scanning model [1,10] and are necessary to project the satellite observation to its correct location on Earth.The exact time, roll, and yaw parameters of the satellite can slightly deviate from these models, causing geometric distortions once the data have been georeferenced based on the provided TLEs.This problem is well known and has been addressed to by several researchers during the decades.Different solutions exist to overcome the imprecise geolocation, both approaching the sensor model itself [11,12] or correcting the already distorted data [13] (more information about these techniques will be presented in Section 3).TIMELINE (TIMe Series Processing of Medium Resolution Earth Observation Data assessing Long-Term Dynamics In our Natural Environment) is a project of the German Remote Sensing Data Center (DFD), German Aerospace Center (DLR).Its aim is to process time series of AVHRR data since the early 1980s.All data must pass a preprocessing, producing a homogenous and intercalibrated basis for the deviation of various level 2 and level 3 products such as vegetation parameters, snow and ice, water bodies, land surface temperature, fire hot spots, and clouds.Time series processing and analysis of these products requires a precise geolocation.Therefore, a method had to be included within the workflow of the TIMELINE processing environment to improve the geolocation accuracy of the data.The proposed two-step chip matching approach is therefore implemented as part of the overall processing flow.In this context, the term "two-step" refers to the back-to-back execution of a matching with image chips extracted from shorelines in the first, and Normalized Difference Vegetation Index (NDVI) composites in the second step, respectively.Figure 1 presents a flowchart of the processors and some dependencies on the product level.It shows that the chip matching is applied prior to all thematic products.This article will summarize the different steps necessary to successfully improve the geolocation accuracy of AVHRR data.Challenges, limitations, and validation results will be presented, offering a detailed insight into the various aspects that need to be considered when planning to implement a processor similar to the one introduced here.

Data
The satellite data processed with the chip matching module originate from all generations of AVHRR (see Table 1).The HRPT data are transferred to the processing chain of the TIMELINE preprocessor from DLRs in-house satellite data archive in raw format.Even though the matching procedure could also be applied to any other sensor, it was designed specifically for the orbit format of AVHRR: the orbit segments that are fed to the procedure can be of variable length (number of scan lines may differ between several hundred and 5600) but the length of each scan line is expected to be 2048 [1].Even though it is designed for HRPT data, Local Area Coverage (LAC) data could also be processed with the presented chip matching approach.The only difference is the record length, which in case of LAC data is limited to ten minutes.
Two reference datasets are required in order to apply the chip matching to an input AVHRR orbit segment: A reference water mask including reference chip positions and as a second reference, a NDVI mask with such chip positions as well.The water mask was created using a dataset provided by naturalearth.com(based on World data bank 2 inventory and updated with a pan-European River and Catchment Database provided by the Joint Research Centre, JRC [14], from 2007) which features a resolution of 10 m and contains all permanent water bodies of the Earth.This dataset was converted to a raster and downscaled to the 0.01 degree resolution of AVHRR.A 20-km buffer zone was added to all coastlines.This is a necessary step for the downstream processing of the water mask that will be used for the matching process.Furthermore, around 400 chip positions were identified along the shorelines of oceans, seas, and large inland lakes.A location qualifies as a chip if the shoreline features a pronounced shape that is going to be unmistakably distinguishable from each possible shift in X-or Y-direction.An example for an unsuitable shoreline would be a linear shape which would look similar if shifted to the north and/or south.The reference water mask is converted to lat/lon WGS84 projection and transformed to orbit projection (the geometry of the original satellite observation from its orbit) each time an AVHRR orbit segment is processed.
The NDVI masks are prepared in a similar way as the reference water mask: The resolution is scaled down to 0.01 degree; and the projection is set to Geographic WGS84.The NDVI information originates from the MODIS operational Vegetation Indices Monthly L3 Global 1 km dataset MOD13 [15].All available monthly products until August 2015 were downloaded and aggregated to create a median NDVI from the 15-year time series of MODIS.Based on this time series, the median NDVI was calculated for each month.The process then scanned quadratic subsets of 64 pixels lateral length which were further split into four 32 × 32 large sub-quadrants, selecting those as chips that feature a significant variance of NDVI between the four sub-quadrants of the subset.For each of the twelve median NDVI datasets, between 1350 and 1650 useful chip locations were detected that way.The number of chips varies as illumination conditions (due to polar night) and the vegetation status differs between seasons.The rules for a useful chip are the same as for the reference water mask: a chip has to be "unique" in its pattern when shifted around, resulting in distinct correlation coefficients when compared to the AVHRR NDVI later in the process.A pattern of a distinct NDVI contrast within each chip is required.The chip locations for the NDVI mask were identified automatically, labeled with a progressional index (each chip gets a unique number) and a reference mask was extracted for each calendar month containing these chips and indexes.

Methods
As indicated in the introduction, there are basically two different approaches for addressing the imprecise geolocation accuracy of AVHRR data: One is to modify the sensor model itself, altering parameters such as observation time, sensor altitude and yaw.This would result in modified TLEs which ultimately would lead to a more accurate geolocation accuracy as demonstrated by, e.g., [11,12].The second possibility is to apply the compromised TLEs that already come with the raw data and then try to correct the resulting pixel shifts with hindsight as it was applied to the Canadian AVHRR Processing System by [16].
Both techniques have their advantages and drawbacks: While the modifications on the sensor model itself result in more precise geolocation accuracy especially under partly clouded conditions, it is difficult to incorporate flight maneuvers executed during the overpass.That is because the necessary modification parameters are extracted from several consecutive observations.The sensor model is adjusted based on the collection of these parameters which originate from a long time span.Orbit trajectory, height, sensor roll and yaw change during a flight maneuver, and it would take a while for the matching process to adjust to the updated conditions.That is not the case for the second approach, which tries to improve the imprecisely georeferenced data and can be conducted without prior knowledge of the sensor model.It attempts to fix the distortions caused by the model rather than to fix the model.The drawback is that it depends on a sufficient amount of precisely identified reference points on the ground and therefore, it requires (at least partially) cloud free observations during daytime.However, as it calculates new matching parameters for each observation, the observations are independent from each other and possible errors do not propagate.
For the TIMELINE project, it was decided to rely on the second method-the subsequent correction of geolocation errors caused by imprecise sensor model parameters following a chip matching approach.The whole geocorrection module is embedded in the preprocessing workflow of the TIMELINE project.It immediately follows the module for calibration, before orthorectification, and is designed as a two-step chip matching approach where the orbit data are compared with reference water and NDVI masks in many different chip positions.For each position, a vector is derived for shifting the orbit data in two-dimensional space.Vectors extracted from the comparison with the reference water mask are independent with those gained from the NDVI comparison.This two-step approach allows for a higher certainty of the resulting matching vectors.The collection of these vectors is finally analyzed to derive the coefficients for a polynomial transformation of the orbit data.Figure 2 contains a schematic workflow of the whole process.The following sections will outline the details of all necessary steps.
Remote Sens. 2017, 9, 303 4 of 17 [11,12].The second possibility is to apply the compromised TLEs that already come with the raw data and then try to correct the resulting pixel shifts with hindsight as it was applied to the Canadian AVHRR Processing System by [16].
Both techniques have their advantages and drawbacks: While the modifications on the sensor model itself result in more precise geolocation accuracy especially under partly clouded conditions, it is difficult to incorporate flight maneuvers executed during the overpass.That is because the necessary modification parameters are extracted from several consecutive observations.The sensor model is adjusted based on the collection of these parameters which originate from a long time span.Orbit trajectory, height, sensor roll and yaw change during a flight maneuver, and it would take a while for the matching process to adjust to the updated conditions.That is not the case for the second approach, which tries to improve the imprecisely georeferenced data and can be conducted without prior knowledge of the sensor model.It attempts to fix the distortions caused by the model rather than to fix the model.The drawback is that it depends on a sufficient amount of precisely identified reference points on the ground and therefore, it requires (at least partially) cloud free observations during daytime.However, as it calculates new matching parameters for each observation, the observations are independent from each other and possible errors do not propagate.
For the TIMELINE project, it was decided to rely on the second method-the subsequent correction of geolocation errors caused by imprecise sensor model parameters following a chip matching approach.The whole geocorrection module is embedded in the preprocessing workflow of the TIMELINE project.It immediately follows the module for calibration, before orthorectification, and is designed as a two-step chip matching approach where the orbit data are compared with reference water and NDVI masks in many different chip positions.For each position, a vector is derived for shifting the orbit data in two-dimensional space.Vectors extracted from the comparison with the reference water mask are independent with those gained from the NDVI comparison.This two-step approach allows for a higher certainty of the resulting matching vectors.The collection of these vectors is finally analyzed to derive the coefficients for a polynomial transformation of the orbit data.Figure 2 contains a schematic workflow of the whole process.The following sections will outline the details of all necessary steps.

Delineation of Cloud Mask, Water Mask, and NDVI from Orbit Data
The matching between reference and orbit data relies on two steps: First, a reference water mask is compared with a water mask derived for an orbit segment to identify a possible shift of the shoreline in as many different chip positions as possible.Second, the NDVI is calculated for the orbit segment and compared with the reference NDVI for the respective month again for all cloud-free chip positions.It was decided to implement this two-step approach for TIMELINE because using two independent methods improves the reliability of the retrieved vectors.
The first step even before a water mask or NDVI can be derived from an observation is the delineation of a cloud mask.At DLR, cloud cover parameters are extracted from AVHRR data relying on the AVHRR Processing scheme Over cLouds, Land and Ocean (APOLLO), an algorithm that was designed in the late 1980s and continuously improved and extended since [2,[17][18][19].Recently, a new version of APOLLO (APOLLO_NG) was developed, now including a probabilistic interpretation of cloud contamination for each pixel within a scene [3].APOLLO provides a huge set of cloud parameters, including optical depth and effective radius, cloud top phase, temperature, and cloud water path.For the chip matching, only the presence of cloud cover is of interest.Therefore, only the basic cloud tests that were already operational in the first APOLLO version (for details see [17]) are implemented in the processing chain and each pixel containing cloud contamination is masked out.
Cloud shadow is another problem that is caused by cloud cover.As the chip matching approach will attempt to match a water mask, and cloud shadow can easily be misclassified as water, cloud shadow has to be dealt with before the water mask is derived.A simple approach was implemented to flag possible cloud cover: Depending on the solar zenith angle (SZA) and the cloud mask, possible cloud shadow is flagged.Cloud height is also an important input variable to the determination of cloud shadow, but this information is not available and therefore, a fixed cloud height (H) of 6 km is assumed.This already includes or excels more than 90% of the anticipated cloud top heights [20].Cloud shadow is approximated under the assumption of a flat earth.The cloud shadow length is then derived following Equation ( 1) and studies published by [21,22].
The direction of the cloud shadow is extracted from the SZA, and the respective area within the orbit data is flagged as possible cloud shadow.Figure 3 illustrates the result of the whole cloud masking approach beginning with the orbit data in a false color composite (band combination 1-2-1 in Figure 3a) and the respective cloud mask (Figure 3b).
The water mask that will be used to match the orbit with a reference water mask is derived from the orbit data relying on a dynamic local threshold determination approach similar to [23,24].The detailed description of the approach is subject of a publication about the TIMELINE Water Mask processor [25].A water mask is derived, using the cloud mask from the first step as input as well as bands 1 and 2 and the reference water mask described in the data section.Within the cloud free pixels of the AVHRR scene, those pixels indicated to be covered by water in the reference water mask are analyzed according to their Near Infrared Reflectance (NIR).A mean NIR is derived from all these pixels as well as the standard deviation (StdDev).For the entire AVHRR scene, the pixels falling within the mean NIR ± 1 StdDev are classified as water, while cloud masked pixels are excluded from this procedure.The procedure has been developed and described by [23].Due to the huge areal coverage of orbital AVHRR data over Europe, the different observation geometry, wave patterns, algae and sediment loads, a fixed threshold classification (e.g., based on the Normalized Difference Water Index [26] or a combination of indices [27]) does not yield satisfying results.Thus, the mean NIR reflectance over cloud free reference regions serves as classification basis.Figure 4 shows the results of the processing steps so far: the calibrated satellite data in false color composite (band combination 1-2-1) on the left (Figure 4a), and the cloud and water masks on the right (Figure 4c).The figures are depicted in the original satellite orbit projection.Finally, the Normalized Difference Vegetation Index (NDVI) is derived from the bands 1 and 2 according to Equation ( 2):  Finally, the Normalized Difference Vegetation Index (NDVI) is derived from the bands 1 and 2 according to Equation ( 2): Finally, the Normalized Difference Vegetation Index (NDVI) is derived from the bands 1 and 2 according to Equation ( 2): With R b1 and R b2 referring to the reflectance in band 1 and 2, respectively.Cloud covered pixels or those prone to cloud shadow are masked out.It is important to mention at this point that during this stage of the processing, atmospheric or BRDF correction have not yet been applied to the data (compare also with Figure 1).Possible effects originating from these omissions cannot be precluded, but can be considered negligible.

Matching Process between Orbit and Reference Data, Transformation
During the chip matching procedure, the derived water mask and NDVI are compared in as many chip positions as possible with the reference datasets introduced in the Data Section, and a possible pixel shift between orbit and reference is extracted for further processing.A two-step chip matching procedure was designed, comparing the orbit data with both a reference water mask as well as a NDVI composite of the respective month.Both steps will return vectors for the subsequent calculation of the warping coefficients, but the determination of these vectors is executed independently for water mask and NDVI.
For each chip within the reference water mask and the NDVI mask, the overall cloud cover percentage is extracted from the respective orbit dataset.If the cloud cover percentage does not exceed 8%, the chip is approved for the matching process and subsequently shifted around for up to 20 pixels in each direction.Figure 5 shows the beginning of the process for a NOAA-17 scene from May 2009.The reference water mask (Figure 5b) is transformed to orbit projection in order to be compared with the orbit water mask (Figure 5c).
With Rb1 and Rb2 referring to the reflectance in band 1 and 2, respectively.Cloud covered pixels or those prone to cloud shadow are masked out.It is important to mention at this point that during this stage of the processing, atmospheric or BRDF correction have not yet been applied to the data (compare also with Figure 1).Possible effects originating from these omissions cannot be precluded, but can be considered negligible.

Matching Process between Orbit and Reference Data, Transformation
During the chip matching procedure, the derived water mask and NDVI are compared in as many chip positions as possible with the reference datasets introduced in the Data Section, and a possible pixel shift between orbit and reference is extracted for further processing.A two-step chip matching procedure was designed, comparing the orbit data with both a reference water mask as well as a NDVI composite of the respective month.Both steps will return vectors for the subsequent calculation of the warping coefficients, but the determination of these vectors is executed independently for water mask and NDVI.
For each chip within the reference water mask and the NDVI mask, the overall cloud cover percentage is extracted from the respective orbit dataset.If the cloud cover percentage does not exceed 8%, the chip is approved for the matching process and subsequently shifted around for up to 20 pixels in each direction.Figure 5 shows the beginning of the process for a NOAA-17 scene from May 2009.The reference water mask (Figure 5b) is transformed to orbit projection in order to be compared with the orbit water mask (Figure 5c).5d-f illustrates the positions of some of the image chips in greater detail.The subset is nearly cloud-free, and the derived orbit water mask (Figure 5f) contains several promising chips (marked with red polygons).The green rectangle from Figure 5f is highlighted in even greater detail in Figure 6.In Figure 6a, the true position of the Ibiza Island is outlined with a green polygon.This information is given within the reference water mask.The position of the island in the orbit data is shown in the background.It is obviously slightly shifted to the south.
Figure 5d-f illustrates the positions of some of the image chips in greater detail.The subset is nearly cloud-free, and the derived orbit water mask (Figure 5f) contains several promising chips (marked with red polygons).The green rectangle from Figure 5f is highlighted in even greater detail in Figure 6.In Figure 6a, the true position of the Ibiza Island is outlined with a green polygon.This information is given within the reference water mask.The position of the island in the orbit data is shown in the background.It is obviously slightly shifted to the south.Figure 6b now only focuses on the area within one chip, where each pixel has been recoded to binary information (1 or 0) about the presence of water.These data are shifted around within the grid as indicated by the blue arrows in Y-and X-direction.Shifts are conducted stepwise, beginning with only one pixel and adding up to twenty pixels in any direction and Y-and X-combination.After each shift the correlation is calculated from the proportion of true (land equal land or water equal water) pixel locations opposed to wrong locations.The combination of X-and Y-shifts resulting in the best correlation is the vector for this chip and will be transferred to an internal database, storing the respective vector as well as the position of the chip within the local grid.This vector is later used for the calculation of the polynomial for the warping process.A vector is only included if the correlation (the linear Pearson correlation coefficient between the binary reference and orbit chip mask) is higher than 0.8 and if the vector is explicitly the only one that results in this correlation (if there is more than one vector result in the same correlation between reference and orbit data, no decision can be made on which vector should be chosen).
The extracted vectors are analyzed to remove outliers by calculating the mean X-and Y-shift within a subset of 1000 × 2048 pixels.These subsets overlap by 200 pixels while moving through the scene, finding and removing the outliers.Additionally, the vectors found for the 150 first and last pixels of a line are removed as the satellite view angle in this region is too wide (approximate ±55.4° to ±50.4° from the nadir [1]) to be able to clearly identify land-water boundaries without blurring effects.If after this step not enough vectors remain for an orbit, the whole process is cancelled.The minimum amount of successfully identified vectors can be calculated using Equation (3) [28]: where n refers to the minimum number of vectors, and p stands for the order of the polynomial warping.In the case of a third order transformation, 10 is the minimum amount of vectors, though it Figure 6b now only focuses on the area within one chip, where each pixel has been recoded to binary information (1 or 0) about the presence of water.These data are shifted around within the grid as indicated by the blue arrows in Y-and X-direction.Shifts are conducted stepwise, beginning with only one pixel and adding up to twenty pixels in any direction and Y-and X-combination.After each shift the correlation is calculated from the proportion of true (land equal land or water equal water) pixel locations opposed to wrong locations.The combination of X-and Y-shifts resulting in the best correlation is the vector for this chip and will be transferred to an internal database, storing the respective vector as well as the position of the chip within the local grid.This vector is later used for the calculation of the polynomial for the warping process.A vector is only included if the correlation (the linear Pearson correlation coefficient between the binary reference and orbit chip mask) is higher than 0.8 and if the vector is explicitly the only one that results in this correlation (if there is more than one vector result in the same correlation between reference and orbit data, no decision can be made on which vector should be chosen).
The extracted vectors are analyzed to remove outliers by calculating the mean X-and Y-shift within a subset of 1000 × 2048 pixels.These subsets overlap by 200 pixels while moving through the scene, finding and removing the outliers.Additionally, the vectors found for the 150 first and last pixels of a line are removed as the satellite view angle in this region is too wide (approximate ±55.4 • to ±50.4 • from the nadir [1]) to be able to clearly identify land-water boundaries without blurring effects.If after this step not enough vectors remain for an orbit, the whole process is cancelled.The minimum amount of successfully identified vectors can be calculated using Equation (3) [28]: where n refers to the minimum number of vectors, and p stands for the order of the polynomial warping.In the case of a third order transformation, 10 is the minimum amount of vectors, though it is recommended to identify at least 15 [28].During the development of the procedure, it was concluded that the minimum number of vectors has to be 18 in order to achieve consistently satisfying results.What follows next is the calculation of additional, artificial vectors within a predefined grid (with an adjustable grid size of 200 pixels).This is necessary as real vectors can only be derived along the shoreline and for some inland lakes (as well as during the NDVI matching, which is an independent procedure explained at the end of this section).Huge areas within the scene will probably not be covered by shorelines, and others might be covered by clouds, adding to the no-data problem for the vector collection.As the intention is to apply a third order polynomial transformation to improve georeferencing at the end, huge areas without any vectors can lead to undesirable distortions-especially if these areas are situated at the corners of the scene.This is why these areas are filled with artificially calculated vectors.To determine the value of an artificial vector, the three spatially closest real vectors are used.Their X-and Y-shift values will be weighted according to the distance to the position of the artificial vector following Equation (4): where WF i is the weighting factor for the real vector i, dist i is the distance of real vector i to the artificial vector position, and n is the number of real vectors being used to calculate the artificial vector, running from 1 to 3. The closer a real vector i is to an artificial vector, the higher its weighting factor WF i .This factor will be determined for all three neighboring vectors, and the final X-and Y-shift of the artificial vector will then be calculated according to Equations ( 5) and ( 6): where Shift x and Shift y are the new vector shifts in X and Y direction, respectively; WF n is the weighting factor for vector n determined from Equation (4); and X n and Y n are the shifts in X-and Y-direction for real vector n, respectively.n is again running from 1 through 3 for the three spatially closest real vectors already used in Equation ( 4).X and Y always refer to the horizontal and vertical position within the AVHRR scene grid, respectively.The above approach results in a collection of vectors (both real and artificial) derived from the shoreline matching.The same procedure is conducted relying on the median NDVI for the respective calendar month as reference data.A similar technique was used by Khlopenkov and Trishchenko [16] to achieve subpixel georeferencing accuracy for AVHRR data in Canada.Even though the reference dataset is completely different from the reference water mask, the methodology is the same: image chips from the reference NDVI are compared with the respective area from the orbit NDVI, the correlation is calculated for each shift in X-and Y-direction, and the collection of real vectors is extended by artificial vectors for regions without sufficient information.
After both matching procedures are finished, the vectors are transferred to the polynomial warping module.Using a least squares estimation, the needed coefficients Kx and Ky of the Equations ( 7) and (8) are determined.
x i = ∑ i,j Kx i,j ×Xo j ×Yo i (7) Ky i,j ×Xo j ×Yo i (8) Xo and Yo contain the original location of the pixels used to determine the coefficients Kx and Ky that are used to map the coordinates to x i and y i .The polynomial warping is performed based on Equation ( 9) and the calculated coefficients: where g[x, y] is the original position of the pixel at coordinates (x,y), while x and y is the pixel that is used to perform the warping process.x and y are calculated by Equations ( 10) and ( 11): The polynomial warping is applied to the latitude and longitude layers of the processed orbit.A quality check is automatically executed comparing the amount of erroneously located pixels along the coastline before and after the matching process.Only if the procedure led to an improvement of the accuracy, the results are actually exported from the workflow and transferred to the subsequent processing chain of TIMELINE.

Validation
To check the validity of the approach, two validation steps were performed: For the first approach, we chose various orbits that originated from different observation dates, comprising all seasons and covering all possible satellite paths over Europe (east, middle and west) and platforms.These orbits were georeferenced using the TLEs provided within the raw data, but no chip matching was applied.Afterwards, we applied the chip matching to the same orbits.The data pairs (unmatched and matched) allow for an evaluation of pixels which were located on the wrong position before the matching and could be improved applying the chip matching.The analysis was performed comparing the distance of the orbit shorelines to the shorelines within the reference water mask.Within a 20 km buffer zone around the shorelines the percentage of false pixels was calculated for matched and unmatched datasets.The direct comparison between these error estimates allows approximating the efficiency of the matching algorithm.
In the second approach, we intentionally modified the latitude and longitude information of the original orbit data artificially to decrease the accuracy of georeferencing.This was done stepwise, starting with only small (one pixel, or 0.01 • ) deteriorations in the first step, and then increasing these artificial errors until a shift of ten pixels was reached.Only scenes with very high initial geolocation accuracy were used for this validation approach.Now the error is known beforehand, and the chip matching can be applied and evaluated afterwards.Not only does this analysis provide an insight into the accuracy of the matching approach, it also offers detailed information about the performance in dependence to the extent of the geometric distortions.The comparison between original and artificially modified data can be performed for randomly distributed pixels and is not limited to the shorelines.The random pixels receive a unique value out of range of the possible measurements and can be traced through the entire chip matching process.The results of both validation steps are presented in the next section.

Results
The described procedure to improve geolocation accuracy relying on a two-step chip matching approach was applied to nearly 600 HRPT datasets.Table 2 contains statistical information about mean, minimum and maximum cloud cover percentage, the standard deviation of the transformation coefficients for X-and Y-direction, and the percentage of daytime pixels (solar zenith angle <85 • ) for both successfully matched orbit segments and those where a chip matching could not be applied.For the statistics presented in Table 2, nearly 600 scenes were processed, 263 of which could be processed successfully.A matching is considered successful if vectors could be extracted, and if the subsequent warping lead to an improvement of the geolocation accuracy.It is possible that no vectors can be identified due to cloud cover or darkness, and it is possible that even though vectors could be extracted, the warping did not lead to an improvement.A reason for such a situation can be that the accuracy already was very high before the matching was attempted.
The comparison between successfully matched and unmatched orbits revealed that in general, cloud cover percentage of the successfully matched orbits is eight percent lower than of the unmatched scenes.In addition, orbits with a higher percentage of nighttime pixels tend to pass through the chip matching without processing.Both effects were expected, as fewer clear-sky pixels automatically lead to fewer unobscured chips.The coefficients that were identified vary basically only in the first position.The standard deviation for the first coefficient is 2.05 for the X-and 1.95 for the Y-direction.Only the first coefficient shows significant variance.The distribution of the single values for the first coefficient is illustrated in Figure 7.For the statistics presented in Table 2, nearly 600 scenes were processed, 263 of which could be processed successfully.A matching is considered successful if vectors could be extracted, and if the subsequent warping lead to an improvement of the geolocation accuracy.It is possible that no vectors can be identified due to cloud cover or darkness, and it is possible that even though vectors could be extracted, the warping did not lead to an improvement.A reason for such a situation can be that the accuracy already was very high before the matching was attempted.
The comparison between successfully matched and unmatched orbits revealed that in general, cloud cover percentage of the successfully matched orbits is eight percent lower than of the unmatched scenes.In addition, orbits with a higher percentage of nighttime pixels tend to pass through the chip matching without processing.Both effects were expected, as fewer clear-sky pixels automatically lead to fewer unobscured chips.The coefficients that were identified vary basically only in the first position.The standard deviation for the first coefficient is 2.05 for the X-and 1.95 for the Y-direction.Only the first coefficient shows significant variance.The distribution of the single values for the first coefficient is illustrated in Figure 7.As expected, no clear pattern is recognizable from the distribution in Figure 7. Shifts can occur in any direction but are concentrated within ±2 in X-and Y-directions, respectively.We also analyzed the resulting vectors over time to see if they are correlated to the date of observation, but they showed no significant dependency in any direction.Only the percentage of successfully processed orbit segments drops during winter, as the solar zenith angles are generally higher and the area percentage affected by polar night is higher, too.

Validation results
The approach to validate the two-step chip matching comprises two parts: The comparison of a reference water mask with AVHRR scenes that have been matched and left unmatched, respectively.For the second part of the validation, geolocation errors have been simulated by artificially shifting As expected, no clear pattern is recognizable from the distribution in Figure 7. Shifts can occur in any direction but are concentrated within ±2 in X-and Y-directions, respectively.We also analyzed the resulting vectors over time to see if they are correlated to the date of observation, but they showed no significant dependency in any direction.Only the percentage of successfully processed orbit segments drops during winter, as the solar zenith angles are generally higher and the area percentage affected by polar night is higher, too.

Validation results
The approach to validate the two-step chip matching comprises two parts: The comparison of a reference water mask with AVHRR scenes that have been matched and left unmatched, respectively.For the second part of the validation, geolocation errors have been simulated by artificially shifting the AVHRR scenes around in X-and Y-directions and then processing these scenes with the chip matching algorithm.Table 3 contains some statistics of the first part: matched vs. unmatched scenes compared with the reference mask.Within the 20 km buffer zone around shorelines, pixels were searched for that were situated in the wrong location (water in the reference mask, land in the orbit data, or vice versa; cloud covered areas have been ignored).In total, approximately 4.18% of all investigated pixels were erroneous, which corresponds to a dislocation of 1.6 km.Applying the chip matching procedure could successfully correct 42.82% of these pixels.The results of the second part of the validation are illustrated in Figure 8, showing the sum of the absolute values of the artificial shift on the X-axis and the percentage of errors introduced by these shifts that could successfully be cleared on the Y-axis.The sum of shifts reaches a maximum of 20 as it is calculated from a maximum X-shift and Y-shift of ten pixels, respectively.Figure 8  the AVHRR scenes around in X-and Y-directions and then processing these scenes with the chip matching algorithm.Table 3 contains some statistics of the first part: matched vs. unmatched scenes compared with the reference mask.Within the 20 km buffer zone around shorelines, pixels were searched for that were situated in the wrong location (water in the reference mask, land in the orbit data, or vice versa; cloud covered areas have been ignored).In total, approximately 4.18% of all investigated pixels were erroneous, which corresponds to a dislocation of 1.6 km.Applying the chip matching procedure could successfully correct 42.82% of these pixels.The results of the second part of the validation are illustrated in Figure 8, showing the sum of the absolute values of the artificial shift on the X-axis and the percentage of errors introduced by these shifts that could successfully be cleared on the Y-axis.The sum of shifts reaches a maximum of 20 as it is calculated from a maximum X-shift and Y-shift of ten pixels, respectively.Figure 8   For the results presented in Figure 8, a total number of 1320 chip matching simulations were performed, 459 of which were successful.Reasons for an unsuccessful cycle include the inability to identify a sufficient amount of unique vectors (~83% of the cases), and the stop of the procedure in For the results presented in Figure 8, a total number of 1320 chip matching simulations were performed, 459 of which were successful.Reasons for an unsuccessful cycle include the inability to identify a sufficient amount of unique vectors (~83% of the cases), and the stop of the procedure in case the matching would lead to a negative result (producing more errors than solving, ~17% of the cases).As explained in Section 3.2, the procedure stops in these events.
The simulation was derived for more scenes, though in these cases not the full range of 440 shifts was completed due to processing time constraints.Thirty-five random combinations of X-and Y-shifts have been applied to 23 additional scenes.The results are different for each case, as was the case for the three scenes illustrated in Figure 8.As indicated above, the chip matching is prone to many external factors such as cloud amount, shape, position, or other data gaps.Therefore, the results of two different scenes will always vary from each other.
The analysis revealed that the chip matching approach presented here is capable of dealing even with large scale geometric errors of up to 10 km, though the effectiveness varies significantly between different scenes.

Discussion
The method to correct the sometimes inaccurate geolocation of AVHRR data relies on a chip matching approach that tries to improve the errors caused by imprecise TLEs.Every observation is processed separately, which causes also independent matching parameters for each scene.The advantage is that errors that can possibly occur when analyzing an observation will not affect yet another scene.The downside of this approach is that, in order for the chip matching to be applied successfully, enough cloud free chip positions must be identified.The presence of cloud cover can vary significantly on a day to day basis, with average cloud cover fractions between 40% (Italy) and 70% (Sweden) [29].Time of the day and the date also influence cloud cover statistics.Therefore, the success rate of the chip matching approach is difficult to predict.A maximum of 8% cloud cover fraction within a chip is allowed, offering some flexibility.Nevertheless, the analysis of the statistics that were generated during processing revealed that around 83% of unsuccessful matching approaches occurred due to cloud cover.For the remaining 17% of unsuccessful approaches, the chip matching did not result in an improvement of the geolocation even though enough chips were identified.
Table 4 Represents a summary of the most important statistics being discussed in this section, giving an overall overview of the capabilities and limitations of the presented chip matching approach.The presented chip matching approach could be applied successfully to 44% of the tested AVHRR scenes (46.5% unsuccessful matchings due to cloud cover, 9.5% of the scenes already provided good geolocation accuracy).This performance might appear amendable on the first view, but scenes with high cloud cover fractions are of little use for time series analysis of land products and, therefore, the inability to improve the geolocation accuracy of these compromised scenes is acceptable.In fact, the analysis of the matching statistics revealed that even within the successfully matched scenes the mean cloud cover percentage was 54% (compare Table 2).This proves that cloud cover itself does not preclude successful matching: it is the position of the clouds that can limit the number of extracted vectors.Additionally, a matching process is only considered successful if the initial geolocation accuracy (the extent of the satellite scene's dislocation with regards to the reference dataset geolocation) could be improved.If this initial accuracy is already high, an additional matching is not necessary (the target geolocation accuracy within the TIMELINE project is one pixel.Thus, an improvement beyond this point is not mandatory).It is difficult to separate these unsuccessful matching cases for several reasons: First, if clouds prevent a successful matching, it is impossible to assess the quality of the initial geolocation.These scenes will be listed under the rubric "unsuccessful due to extensive cloud cover" even though it is possible that the quality is already good enough anyway.Secondly, it is very unlikely that all shifting vectors that are extracted from a cloud-sparse scene will hold values of only zeros in X-and Y-direction.The distortions in the off-nadir region of an observation tend to induce at least small shifts which result in subpixel adjustments of the geolocation.It is important to emphasize that also subpixel adjustments to improve geolocation accuracy are feasible with the two-step chip matching approach presented here: As the final warping is applied to the values within the lat/lon layer even modifications of few meters only can be realized.This effect is illustrated in Figure 7 by the concentration of X-and Y-coefficients around-but not exactly-zero.As a consequence of one or more vectors uneven zero, a polynomial transformation is applied to the lat/lon layer, and the internal quality control will assess whether or not the geolocation accuracy could be improved.In case of only very slight subpixel distortions in the off-nadir region, the transformation does not always lead to an improvement of the geolocation accuracy of the whole scene, and the results are therefore discarded.As described in Section 4, such incidents occur in 17% of unsuccessful matching approaches.
The accuracy assessment was derived following two approaches: In the first step pairs of water masks extracted from unmatched and matched AVHRR scenes were compared with a reference water mask to identify the amount of erroneous pixels.This investigation revealed that approximately 4.2% of all pixels within a 20 km buffer around the shoreline originating from the unmatched scenes were incorrect.This would result in an average geolocation uncertainty of around 1.6 km.Comparing the matched water masks with the reference water mask showed that 42.8% of the existing errors could be corrected-for those scenes that could successfully be processed by the two-step chip matching.This results in an average remaining error of 0.92 km, which is below the target geolocation accuracy of one pixel (~1.1 km).Again, the improvement of 42.8% does not seem extraordinarily high on a first look (as it was the case with the success rate of 44%), and again the main reason for this outcome is cloud cover.The more cloud-free chips allow for vector extractions, the more inputs can be provided to the polynomial transformation.The average improvement of 42.8% includes scenes with cloud cover percentages ranging from 22% to 81% (compare Table 2).Understandably, the polynomial extracted for a scene with 81% cloud cover will not lead to an improvement as high as for a scene with only 22% cloud cover.
The second validation approach was performed simulating pixel shifts in several AVHRR scenes and trying to repair these artificial geolocation errors afterwards.This process turned out to be computationally intensive: We simulated shifts of up to ten pixels in any direction and any combination, which lead to 440 processes for each scene (−10 through 10 in both X-and Y-directions, while the combination 0-0 was not processed as it represents the reference for the comparison).Figure 8 illustrates the percentage of errors that could be repaired depending on the sum of the absolute values of the artificial pixel shifts for three example scenes.The graphs show two characteristic properties of the presented chip matching approach: The percentage of improved pixels is lowest for scenes with only slight geometric distortions.According to the second validation approach, the effectiveness of the matching process increases with growing geometric imprecision: While the effectiveness reaches only between 10% and 30% percent for shifts between −1 and 1 in any direction (X and/or Y), it can increase to more than 60% for pixel shift sums of more than 15.This leads to the conclusion that larger deviations can be compensated more effectively than very small geometric disparities.The reason is that the polynomial transformation is performing a warping of the lat/lon layer according to vectors which are valid only for the respective image chips.The area between these chips is interpolated, which introduces small inconsistencies.If the initial overall error is already huge, these inconsistencies are negligible-in relation to the initial error-and the corresponding improvement rate is high.If on the other hand the initial error is very small (e.g., one pixel), the inconsistencies that are introduced in the region between the image chips reduce the overall effectiveness.The second characteristic that can be identified from Figure 8 is that the effectiveness of the matching varies significantly between the analyzed scenes.This is due to the different coverage of the individual orbit segments, the varying amount and location of image chips available for the respective scene, and-as was already mentioned several times-the variable cloud cover conditions.Additionally, the matching with NDVI chips (second step of the two-step approach) may lead to unsatisfactory results if the vegetation in a particular scene deviates considerably from the long term mean for the respective month.This can happen for scenes at the very beginning/end of spring or autumn months, as vegetation activity can differ from long term means at these dates.The algorithm will try to compare the actual AVHRR NDVI with long term mean NDVI derived from MODIS between 2000 and 2015.Various conditions can influence the actual NDVI, such as exceptionally long or short snow cover duration during spring, early or late snow cover onset during autumn, and intense or little rain fall (and therefore onset of vegetation growth).Even though general vegetation patterns will still be present within the designated chips (e.g., forest lines, vegetation/water body boundaries, tree lines within mountainous areas), the overall amount of valid chips might decrease, which affects the effectiveness of the matching.It is also possible that the conditions in some parts of Europe meet the average expectations while in other regions, vegetation growth might run late.All these possible variations influence the performance of the two-step chip matching approach.An internal quality control ensures that the matching is only applied to the lat/lon layers if it actually improves the accuracy.As the subsequent polynomial transformation is then only applied to the lat/lon layer, the original satellite data remeain untouched.No interpolation between pixel values is necessary at this point of processing, preserving the pure physical measurements.Additionally, the polynomial that will be used to improve geolocation accuracy will be incorporated within the metadata of the final Level 1B-product.A quality layer is provided which informs about the reliability of every 512 × 512 pixel subregion of an image.Potential users can analyze this layer, e.g., if they are only interested in products with a high certainty that geolocation is accurate for a specific area within Europe.
The two-step chip matching approach is designed in a way that it can be transferred to other regions outside Europe.A reference water mask containing image chips as well as NDVI layers for each month are required, but they can be produced easily from any arbitrary source.The orbit length (number of lines) is flexible, which allows processing in any desired region.Only if sensors other than AVHRR are intended to be processed by the presented algorithm, slightly extended modifications would become necessary as the number of rows is a fixed constant for AVHRR, and some classification thresholds would require revisions.However, since the matching is performed relying on water masks and NDVI, the method can generally be applied to any desired sensor capable of extracting water masks and calculating NDVI from its measurements.

Conclusions
The two-step chip matching algorithm presented here was designed to improve the geolocation accuracy of Advanced Very High Resolution Radiometer (AVHRR) data, which are the basis for the TIMELINE (TIMe Series Processing of Medium Resolution Earth Observation Data assessing Long-Term Dynamics In our Natural Environment) project at the German Remote Sensing Data Center (DFD), German Aerospace Center (DLR).The target geolocation accuracy within TIMELINE is one pixel (~1.1 km), which is not always provided by the original AVHRR data.
The matching is performed fitting subsets (image chips) of a reference water mask with a water mask extracted from the satellite data in a first step.For each successfully fitted chip, a vector is extracted indicating the shift that is necessary to improve the geolocation accuracy.In the second step, image chips originating from a long-term median NDVI reference mask are fitted to the NDVI calculated from the satellite data.Again, shifting vectors are extracted for each successful match.The vectors from both matching steps are analyzed for consistency, gaps between the vector positions are interpolated, and a third order polynomial transformation is applied to the lat/lon layer of the satellite data.The data layers of the AVHRR imagery remain unchanged to preserve the original satellite measurements.An internal quality control ensures that the matching is only applied if the initial geolocation accuracy can be improved.
The accuracy assessment revealed that the average geolocation error within the original AVHRR data is around 1.6 km.This error could be reduced to an average of 0.92 km for all successfully processed scenes.The effectivity of the matching is low for very low geolocation errors, as the extracted vectors are only valid for the respective image chip position.The area between these chips is interpolated which introduces minor uncertainties.If the initial accuracy is already high, the overall improvement is lower: Around 20% of shifts within ±1 pixel could be corrected.This rate can increase to more than 60% if the initial geolocation error is higher.
Cloud cover turned out to be the biggest obstacle.A total of 46.5% of all processed scenes could not be matched because of clouds covering too many of the image chip positions.As cloud cover is very dynamic and varies considerably between different observations, each processed scene differs from its precursors.A total of 9.5% of all processed scenes turned out to feature a comparably high initial geolocation accuracy which could not be improved any further.This results in an overall success rate of the presented matching approach of 44%.
The algorithm can generally be applied to any desired region also outside Europe, as long as reference water and NDVI masks are provided.If sensors other than AVHRR have to be processed, some modifications will be necessary as classification thresholds as well as additional parameters such as number of data rows or pixel size would change.

Figure 1 .
Figure 1.Flowchart of TIMELINE products and dependencies with the position of the chip matching processor highlighted in dark grey.Apollo NG: Next Generation Apollo [3] (cloud processing scheme); BRDF: Bidirectional Reflectance Distribution Function; LST: Land Surface Temperature; NDVI: Normalized Difference Vegetation Index; LAI: Leaf Area Index; FVC: Fractional Vegetation Cover; FAPAR: Fraction of Absorbed Photosynthetically Active Radiation.

Figure 1 .
Figure 1.Flowchart of TIMELINE products and dependencies with the position of the chip matching processor highlighted in dark grey.Apollo NG: Next Generation Apollo [3] (cloud processing scheme); BRDF: Bidirectional Reflectance Distribution Function; LST: Land Surface Temperature; NDVI: Normalized Difference Vegetation Index; LAI: Leaf Area Index; FVC: Fractional Vegetation Cover; FAPAR: Fraction of Absorbed Photosynthetically Active Radiation.

Figure 2 .
Figure 2. Flowchart of chip matching procedure; chip masks include both, the reference water mask as well as an median NDVI mask for the respective month; APOLLO: AVHRR Processing scheme Over Land cLoud and Ocean.

Figure 2 .
Figure 2. Flowchart of chip matching procedure; chip masks include both, the reference water mask as well as an median NDVI mask for the respective month; APOLLO: AVHRR Processing scheme Over Land cLoud and Ocean.

Figure 3 .
Figure 3. Cloud and cloud shadow mask for a NOAA-17 scene from 17 August 2009, depicting: (a) the false color orbit data in band combination 1-2-1; and (b) the cloud mask including cloud shadows in red.The figure is depicted in orbit projection.

Figure 4 .
Figure 4. NOAA-17 Scene of 17 August: (a) in false color band combination 1-2-1; and (b) the cloud and water mask; both figures in orbit projection.

Figure 3 . 17 Figure 3 .
Figure 3. Cloud and cloud shadow mask for a NOAA-17 scene from 17 August 2009, depicting: (a) the false color orbit data in band combination 1-2-1; and (b) the cloud mask including cloud shadows in red.The figure is depicted in orbit projection.

Figure 4 .
Figure 4. NOAA-17 Scene of 17 August: (a) in false color band combination 1-2-1; and (b) the cloud and water mask; both figures in orbit projection.

Figure 4 .
Figure 4. NOAA-17 Scene of 17 August: (a) in false color band combination 1-2-1; and (b) the cloud and water mask; both figures in orbit projection.

Figure 5 .
Figure 5. Example for chip matching procedure: (a) reflection in band 2 (NIR); (b) reference water mask; and (c) derived orbit water mask.The subsets below contain the same features in detail (d,e), with (f) indicating the position of the image chips for the future matching process in red polygons.The green rectangle in (f) marks the area illustrated in even greater detail in Figure 6 (Ibiza Island).All figures in orbit projection

Figure 5 .
Figure 5. Example for chip matching procedure: (a) reflection in band 2 (NIR); (b) reference water mask; and (c) derived orbit water mask.The subsets below contain the same features in detail (d,e), with (f) indicating the position of the image chips for the future matching process in red polygons.The green rectangle in (f) marks the area illustrated in even greater detail in Figure 6 (Ibiza Island).All figures in orbit projection

Figure
Figure5d-f illustrates the positions of some of the image chips in greater detail.The subset is nearly cloud-free, and the derived orbit water mask (Figure5f) contains several promising chips (marked with red polygons).The green rectangle from Figure5fis highlighted in even greater detail in Figure6.In Figure6a, the true position of the Ibiza Island is outlined with a green polygon.This information is given within the reference water mask.The position of the island in the orbit data is shown in the background.It is obviously slightly shifted to the south.

Figure 6 .
Figure 6.Detail from Figure 5f, showing: (a) The position of one single image chip (in red) superimposed on the orbit water mask.The green polygon depicts the shape of the true coastline (as originating from the reference water mask), the cyan polygon is an island outside the actual chip (Formentera Island).(b) The chip area only is illustrated, including a pixel grid and the orbit water mask as a binary (black/white) layer in the background.The blue arrows indicate the four directions in which the binary orbit is shifted around to find the best match.

Figure 6 .
Figure 6.Detail from Figure 5f, showing: (a) The position of one single image chip (in red) superimposed on the orbit water mask.The green polygon depicts the shape of the true coastline (as originating from the reference water mask), the cyan polygon is an island outside the actual chip (Formentera Island).(b) The chip area only is illustrated, including a pixel grid and the orbit water mask as a binary (black/white) layer in the background.The blue arrows indicate the four directions in which the binary orbit is shifted around to find the best match.

Figure 7 .
Figure 7. Distribution of first coefficient values for polynomial warping.The crosses show the X-and Y-shift of all successfully processed orbits.

Figure 7 .
Figure 7. Distribution of first coefficient values for polynomial warping.The crosses show the X-and Y-shift of all successfully processed orbits.
contains three graphs, originating from the simulation of pixel shifts for three different scenes, a NOAA-17 scene from August 2009 (grey graph), a NOAA-19 scene from May 2009 (red graph), and a NOAA-16 scene from June 2001 (blue graph).Remote Sens. 2017, 9, 303 12 of 17 contains three graphs, originating from the simulation of pixel shifts for three different scenes, a NOAA-17 scene from August 2009 (grey graph), a NOAA-19 scene from May 2009 (red graph), and a NOAA-16 scene from June 2001 (blue graph).

Figure 8 .
Figure 8. Ability of chip matching procedure to correct artificially introduced geometric errors for three different scenes.Grey graph is of August 2009, red graph is of May 2009, blue graph is of June 2001.The X-axis represents the sum of the absolute values (modulus) of X-and Y-shifts (artificial errors) within each AVHRR scene, the Y-axis illustrates the percentage of errors introduced by these shifts that could be cleared successfully.

Figure 8 .
Figure 8. Ability of chip matching procedure to correct artificially introduced geometric errors for three different scenes.Grey graph is of August 2009, red graph is of May 2009, blue graph is of June 2001.The X-axis represents the sum of the absolute values (modulus) of X-and Y-shifts (artificial errors) within each AVHRR scene, the Y-axis illustrates the percentage of errors introduced by these shifts that could be cleared successfully.

Table 1 .
Spectral characteristics of Advanced Very High Resolution Radiometer sensors.

Table 1 .
Spectral characteristics of Advanced Very High Resolution Radiometer sensors.

Table 2 .
Statistics of orbit segments processed by chip matching.

Table 2 .
Statistics of orbit segments processed by chip matching.

Table 3 .
Percentage of wrong pixel locations within the buffer zone before the chip matching, and percentage of achieved improvement based on the number of wrong pixel location after applying the matching.

Table 3 .
Percentage of wrong pixel locations within the buffer zone before the chip matching, and percentage of achieved improvement based on the number of wrong pixel location after applying the matching.

Table 4 .
Overview of chip matching statistics.