A Parallax Shift Effect Correction Based on Cloud Height for Geostationary Satellites and Radar Observations

The effect of cloud parallax shift occurs in satellite imaging, particularly for high angles of satellite observations. This study demonstrates new methods of parallax effect correction for clouds observed by geostationary satellites. The analytical method that could be found in literature, namely the Vicente et al./Koenig method, is presented at the beginning. It approximates a cloud position using an ellipsoid with semi-axes increased by the cloud height. The error values of this method reach up to 50 meters. The second method, which is proposed by the author, is an augmented version of the Vicente et al./Koenig approach. With this augmentation, the error can be reduced to centimeters. The third method, also proposed by the author, incorporates geodetic coordinates. It is described as a set of equations that are solved with the numerical method, and its error can be driven to near zero by adjusting the count of iterations. A sample numerical solution procedure with application of the Newton method is presented. Also, a simulation experiment that evaluates the proposed methods is described in the paper. The results of an experiment are described and contrasted with current technology. Currently, operating geostationary Earth Observation (EO) satellite resolutions vary from 0.5 km up to 8 km. The pixel sizes of these satellites are much greater than for maximal error of the least precise method presented in this paper. Therefore, the chosen method will be important when the resolution of geostationary EO satellites reaches 50 m. To validate the parallax correction, procedure data from on-ground radars and the Meteosat Second Generation (MSG) satellite, which describes stormy events, was compared before and after correction. Comparison was performed by correlating the logarithm of the cloud optical thickness (COT) with radar reflectance in dBZ (radar reflectance – Z in logarithmic form).


Introduction
The precision of remote space observations is important when investigating and monitoring various components of global ecological systems, such as marine, forestry, and climate environments [1][2][3][4]. Satellite data integration with external marine and other datasets is crucial in various applications of remote sensing techniques [5,6]. For climate and meteorological investigations, observations of clouds and precipitation on a global scale are usually performed using ground-based radar data and observations from geostationary satellites, due to their high temporal and moderate spatial resolution [7][8][9]. However, during data comparison and integration from these sources, the problem of parallax shift occurs [7,10], which is particularly observable for mid and high latitudes, and also for longitudes far from the sub-satellite point. Parallax shift is also important for cloud shadow determination, which is a significant issue for solar farms [11] and for flood detection [12]. Parallax phenomena also have a significant impact on the comparison of data from low-orbit satellites from different sensors [13][14][15][16].
In terms of mathematical problem formulation, the parallax shift effect for the geostationary satellites is actually a special case amongst low-orbit satellites, and it is easier to investigate due to higher temporal resolution data acquisition and the fixed satellite position.
There have been several attempts to solve parallax shift for geostationary satellites. One of them was proposed by Roebeling et al. [7,17,18], and was based on liquid water path (LWP) value pattern matching. This approach was suitable for stormy events and other inhomogeneous cloud formations, however it usually failed to perform correction in the case of homogeneous spatial LWP distribution. Another attempt proposed by Greuell et al. and Roebeling [19,20] used a simplified geometric model, which assumes Earth to be locally flat, as well as a sort of a priori knowledge about cloud height above the Earth's surface. There were also attempts by Li, Sun, and Yu [12] to solve this problem using a spherical model. Finally, there is the Vicente et al./Koenig method [21,22] based on the geometric properties of parallax shift phenomena, incorporating an ellipsoid model of Earth. The Vicente et al./Koenig method will be presented further in this paper.
There are two methods proposed by the author in this paper, which are based on the same assumptions as the Vicente et al./Koenig method. The first is an augmented version of the mentioned method. This augmentation reduces the correction error to centimeters. The second method proposed by the author is an original work which is based as before on a priori knowledge of cloud height, a geodetic equation of an ellipsoid, and numerical methods for solving the equation set. This method allows the correction error to be reduced to almost zero (assuming Earth to have an ellipsoidal shape).

Problem Description
A parallax shift error in satellite observations occurs when the apparent image of the object is placed in the wrong location on the ellipsoid, considering the ellipsoid's normal line passing through the observed point. This geometric phenomena is particularly observable in geostationary and polar satellite observations due to the high angles of observations, particularly for edge areas of image scenes. In Figure 1, the problem is presented considering the case of a geostationary satellite. As a result, this phenomena causes pixel drift from the original position towards the edge of the observation disk. Consequently, the higher the cloud top layer is, the bigger the shift that occurs.  The position of the cloud top (T) in Cartesian coordinates can be formulated as follows [23]: is the prime vertical radius of the curvature, 2 = 2 − 2 2 is the square of eccentricity, is Earth's semi-major axis, is Earth's semi-minor axis, ℎ is the cloud top height, , is the geodetic latitude and longitude, and 0 is the longitude above which the geostationary satellite is floating. In this case, Equation (1) models the cloud position on a tangent line at coordinates , (see Figure 2). This is a more precise model than the flat-earth model or the spherical model. Note that: all longitudes ( , and ) are equal and the same. Subscripts are given to formally distinguish these values between corresponding latitudes that have different definitions (see Figure  2). is the parametric latitude, P is the point of interest on the ellipsoid, and P* is the image of the point of interest on a sphere. Based on figures from [23,24].
Pixel displacement in satellite view coordinates is defined as: where and are constants that allow for sensor inclination angles to be converted to pixels or distance units in the satellite view space. Also, (ℎ) and (ℎ) are defined as: Downloaded from mostwiedzy.pl where (ℎ), (ℎ), and (ℎ) are cloud top coordinates from Equation (1) as functions of h; = + ℎ -distance from center of Earth to satellite; is the Earth's semi-major axis; ℎ is the distance from the surface of Earth to the satellite. In order to illustrate (ℎ), the following analysis presented in Figure 3 was performed. Namely, depending on the geographical localization of the affected pixel and cloud top height, the absolute shift error in observations is expressed in Spinning Enhanced Visible Infra-Red Imager (SEVIRI) pixel units (In this case = = ℎ 3 ⁄ ). It is worth noting that in many cases, especially for observations of clouds over 5000 m, this can cause pixel shift in the SEVIRI instruments used for the purpose of this study. As mentioned earlier, this effect hinders the comparison process between satellite and groundbased radar data [7]. An example is depicted in Figure 4.

Vicente et al./Koenig Method
The parallax shift problem is solved using a geometrical model, assuming that the surface of Earth is an ellipsoid, and with a priori knowledge of cloud top height. One of the approaches considered in this work is the method proposed by Vicente et al. [21] and implemented by Marianne Koenig [22]. This approach, similar to the rest of the methods presented in this paper, assumes a priori knowledge of the cloud top height, which can be calculated using the observed brightness temperature [7,25]. In this method, the Cartesian coordinates of the cloud image are described as: where is Earth's semi-major axis; is Earth's semi-minor axis; ℎ is the cloud top height; and are the geocentric latitude and longitude (see Figure 2), respectively; 0 is the latitude of the geostationary satellite position; and ( ) is the local radius of ellipsoid for the geocentric latitude model: where: Satellite position (S) is defined as: where is Earth's semi-major axis and ℎ is the distance from the surface of Earth to the satellite.
▪ Radar data ▪ Satellite data Downloaded from mostwiedzy.pl 4. Designate coefficient , which allows Cartesian coordinates of the cloud top to be calculated using the following equations (see Figure 5): where | | ��������⃗ is described by the ellipsoid parametric equation: where and are the parametric latitude and longitude (see Figure 2). Therefore, Equation (9) can be presented as a set of equations: Squaring each equation and adding them according to their sides leads to a square equation, which can be solved with respect to c: 5. Apply c to calculate the Cartesian coordinates of -| | ���������⃗ , | | ���������⃗ , and | | ���������⃗ . 6. Calculate the geocentric ellipsoidal coordinates of : 7. If required for further computation, a geodetic latitude can be calculated: Note that Equation (10) does not describe the cloud top position as it was defined in Equation (1) in Section 2.1. The coordinates of the point are shifted to height ℎ above the ellipsoid along the normal vector. Instead, it describes the point on the ellipsoid with the semi-axes increased by ℎ, therefore this method is burdened with error because of the inadequacy of the model.

Vicente et al./Koenig Augmentation
The Vicente et al./Koenig method can be augmented in the final steps, where the latitude of the cloud bottom position is calculated. When using the Vicente et al./Koenig method, it is assumed that the cloud top is located on the ellipsoid with semi-axes increased by h, and therefore the geodetic latitude can be calculated taking into account the mentioned assumption: If further computation requires the geocentric latitude, this can be calculated using the following equation: This modification allows the correction error to be reduced to centimeters. Details will be presented in the experimental section.

Ellipsoid Model with Geodetic Coordinates: Numeric Method
This method incorporates the cloud top position defined in Section 2.1 in Equation (1). With the described cloud top position, the geostationary satellite observation line should be defined as: where = + ℎ is the distance from Earth's center to the satellite; is Earth's semi-major axis; ℎ is the distance from the surface of Earth to the satellite; and are satellite inclination angles; is the distance from the satellite along the observation line. To solve this problem, an intersection point between the surface above the ellipsoid and the observation line needs to be calculated. Equations (1) and (17)  The root of the above system of equations ( and ) is the geodetic coordinates of point B. However, due to the entanglement of the variable in Equation (18), the root of the equations was designated using the numerical approach. The above method was implemented in C++ and Matlab. The Matlab implementation uses the fsolve function [26], which is part of the optimization toolbox. A detailed configuration of the fsolve function will be presented in the next section. The C++ implementation incorporates the Newton method, which is described below.
To solve the problem using the Newton method, the target function should be defined: where: In Equation (19), � ( , , )� can be interpreted as the distance between the current solution and the optimal solution, which in the optimal case should be equal to zero. For such a defined cost function, calculation of the next iteration of the solution for the Newton method is defined as: where: and: and: The stopping condition is defined as: However, the convergence of the above-presented approach is difficult to obtain for areas located near the edges of the observation disk. Therefore, an alternative target function is defined as the element-wise square of Equation (19): with the gradient defined as: and the stop condition:

Downloaded from mostwiedzy.pl
Another issue that occurs in the numerical calculation problem is the big difference in scale between , , which are expressed in radians, and , which is expressed in meters. To handle this problem, all distances ( , , ℎ, , ) should be divided by . This operation will bring to a similar scale as , . An example result of parallax correction using the numerical method via the Newton algorithm is presented in Figure 6. Note that the radar data is better aligned with the satellite data than in Figure  4. Figure 6. Comparison of detected precipitation mask based on ground-based radar data (blue) and data from Meteosat Second Generation with applied parallax correction using a numerical algorithm (green). Images of small storm clouds from satellite and meteorological radars in the bottom-right corner seem to overlap. The height of the cloud tops reaches 12 km. The stormy event is dated July 25, 2015, 13:00 UTC. EuroGeographics was used for the administrative boundaries.

Parallax Effect Correction Error Simulation
In order to compare the parallax effect correction obtained by the analyzed methods, a simulation experiment was performed. The main goal of the experiment was to generate several cloud top positions that simulate geostationary satellite observations, which result in and for simulated cloud top heights. With the and coordinates and a priori knowledge of the cloud height, correction methods were performed and their results were compared with the original (simulated) cloud position. The detailed procedure of the experiment is as follows: 1. Prepare a grid of geodetic coordinates: ∈ 〈−90°; 90°〉, ∈ 〈−90°; 90°〉, with 1° steps for each dimension; 2. Transform the grid coordinates to the geostationary view coordinates system, , (from now on called the base , ) [27], and back to geodetic coordinates to specify which grid elements are out of scope; for out-of-scope elements, this operation will return Not a Number (NaNfloating point special value Downloaded from mostwiedzy.pl e. The distance between the simulated original base , and ′ , ′ in the geostationary view space will be denoted as the correction error. The correction error is calculated in the geostationary view coordinate space (violet surface on Figure 1), because it allows the impact of the correction error on a specific satellite sensor to be estimated. The coordinates in the above equation were expressed as an angle, however expressing them in radians and multiplying by ℎ allows the result to be calculated in metric units (meters) as distances on a sphere of radius ℎ around a geostationary satellite. This interpretation of geostationary coordinates is implemented in the PROJ software library [28].
In order to calculate the correction using the geodetic coordinates numerical method, the fsolve [26] function was applied. All distances were normalized with respect to the radius of the equator. The parameters of the fsolve function were as follows: • Algorithm: Levenberg-Marquardt (instead of Newton); • Function tolerance: 200 ⁄ ; • Specify objective gradient: yes; • Input damping: 10 −5 . The results of the simulation using the Vicente et al./Koenig method and its augmented version are presented in Figure 7 and 8. The results using the geodetic coordinates numerical method are presented in Figure 9 and 10.
In Figure 7, the errors of the Vicente et al./Koenig method and its augmented version are depicted for certain cloud heights. Note that the error for the augmented version is 10 3 times smaller than for the unmodified version. Also, the median of error rises near linearly with the increase of the cloud height. Also note that the error rises as the distance from the equator and from the central meridian increases. Figure 8 shows histograms of the errors presented in Figure 7. In the histograms, the error ratio between Vicente et al./Koenig and its augmented version can also be spotted, which can be estimated as 10 3 . Another important piece of information is that for the assumed cloud heights, the maximal error of the Vicente et al./Koenig method can be estimated at 50 m, and for the augmented version, it can be estimated at 5 cm.
The errors of the geodetic coordinates numerical method for chosen cloud heights along with the number of iterations of the numerical method are shown in Figure 9. Note that the error is below 1 cm for almost the entire disc. The biggest errors appear near the edges in regions where the Vicente et al./Koenig method failed to compute a result (red NaN regions in Figure 7). The number of iterations increases as the height of the clouds and the distance from the center of the observation disc increase. However, during the performed experiments, the value for the majority of cases was less than or equal to five.
The histograms of errors for the geodetic coordinates numerical method and its number of iterations are presented in Figure 10. Based on the obtained results, the error histograms seem to be quite similar between the experiments-almost all values are classified as near zero. However, there are several occurrences of errors up to 3 meters, which are mainly caused by pixels in regions near the edge of the observation disc. The iteration histograms evolve along with the cloud height. As can be seen, the majority of occurrences fall below five iterations. Occurrences above this value refer to regions near the edge of the observation disc. Downloaded from mostwiedzy.pl

i) j)
Downloaded from mostwiedzy.pl

i) j)
Downloaded from mostwiedzy.pl

i) j)
Downloaded from mostwiedzy.pl

i) j)
Downloaded from mostwiedzy.pl

Discussion
The results of the conducted experiments indicate that the Vicente et al./Koenig parallax effect correction method error is smaller than 50 meters in the geostationary satellite coordinates system for cloud heights of up to 16 km. The error for the augmented version of the Vicente et al./Koenig method proposed by the author allows the error values to be decreased to below 10 cm, which is negligible for current practical applications. As expected, the error of the geodetic coordinates numerical method is also negligible because it can be adjusted by the number of iterations. However, the advantage of the numerical approach is that it corrects the positions of pixels located near the edge of the observation disc (there is no NaN on Figure 9).
On the other hand, it must be noted that the proposed approach requires greater computational power than a method with a constant number of steps, such as the Vicente et al./Koenig method. However, experiments show that this is negligible, as the parallax correction problem was computed within minutes.
As was mentioned in the introduction, parallax effect correction is significant for the comparison and collocation of meteorological radar data and geostationary satellite data. This can be demonstrated by comparing radar reflectance in dBZ: where Z is radar reflectance. Reflectance is described as the following empirical relation with precipitation rate R [mm/h] [29]: and cloud optical thickness (COT) in logarithmic form [30,31]. Figure 11. Scatterplot representing the dependence of radar reflectance and cloud optical thickness (logarithm) data for a stormy event on July 25, 2015, 13:00 UTC, without parallax effect correction (see Figure 4). The calculated Pearson's correlation coefficient is 0.556. Figure 11 presents a scatterplot for radar reflectance and cloud optical thickness for satellite data without parallax effect correction, which should be mutually correlated in ideal cases. The Pearson's Downloaded from mostwiedzy.pl correlation value for that case is equal to 0.556. On the other hand, Figure 12 presents the same type of scatterplot but for the satellite data after parallax effect correction (by numerical method from Section 3.2). The correlation value for the corrected data is equal to 0.683. Note that a threshold effect occurs on top of both figures (presented as a horizontal set of points equal to 2.4), which is a consequence of Optimal Cloud Analysis (OCA) algorithm look-up table (LUT) limitations [31]. It is worth noticing that this effect is less significant for Figure 12, suggesting that data with parallax effect corrections are improved in terms of geometric accuracy. Figure 12. Scatterplot representing the dependence of radar reflectance and cloud optical thickness (logarithm) for a stormy event on July 25, 2015, 13:00 UTC, with parallax effect correction (see Figure  6). The calculated Pearson's correlation coefficient is 0.683.
Note that despite performed spatial correction, data presented in Figure 6 and Figure 12 still differ. In this context, it is important to note that these differences are caused by other factors that influence data acquisition, namely: • Different nature of the acquisition model, as on-ground radar and MSG satellite acquisition are registered with a slight temporal shift (less than 15 min); • Both sensors utilize the different physical natures of acquisition. The on-ground radar is an active sensor which sends out an electromagnetic signal in the microwave spectrum and measurers the echo intensity scattered from precipitation particles. On the contrary, MSG SEVIRI is a passive sensor that measures radiation in a particular electromagnetic bandwidth (visible and near visible spectrum) coming from the sun and thermal radiance; • Data acquired by MSG and on-ground radar is also characterized by different spatial resolutions. Therefore, in order to compare these datasets, additional resampling needs to be performed. Another aspect worth considering is algorithm sensitivity to the uncertainty of cloud top height. The easiest way to approximate this is to calculate the sensitivity of the parallax error itself due to changes of cloud top height. The sensitivity of the parallax error in satellite coordinates is defined as a derivative of pixel displacement (Equation (2)) in respect to h: Downloaded from mostwiedzy.pl Because pixel displacement is nearly linear in respect to h (as can be noticed on Figure 3), the derivative (Equation (31)) should nearly be constant for assumed and . Therefore, it can be approximated as the mean slope (ℎ) in respect to h, for instance: where: = = ℎ . The displacement sensitivity depends on and , therefore its value varies around the globe. Sensitivity values for cities from Figure 3 are presented in Table 1. Note that, the displacement sensitivity can be roughly approximated as less than 1. Therefore, SEVIRI instrument uncertainty of cloud height greater than or equal to 3 km may lead to one pixel size or greater error.

Conclusions
Data integration with data acquired from different sources requires developing additional methods that aim to reduce the discrepancies resulting from different physical aspects of observation. In this context, parallax shift correction for satellite data is a process that reduces geometric differences between observations, and in many cases can significantly improve the quality of corrected data in comparison with on-ground sources.
Regarding the scope of practical applications of the proposed approaches, it is important to note that the resolution of currently operating geostationary satellites varies between 1-3 km for a SEVIRI instrument [32] and 1-8 km for a Geostationary Operational Environmental Satellite (GOES) imager [33]. The upcoming series of Meteosat Third Generation (MTG) satellites will provide imagery with a spatial resolution between 0.5 and 2 km [34]. All of the above-presented parallax methods are effective enough for current and near future geostationary observation satellites, and the usefulness of the proposed methods is negligible. Selection of the proposed parallax effect correction method will be significant only when the spatial resolution of geostationary observations is comparable to 50 Downloaded from mostwiedzy.pl m. With high data resolution and precise parallax effect correction, the algorithm influence of precision on cloud height will become noticeable.
The parallax shift phenomena also affect the comparison between data collected from geostationary satellites and low-orbit satellites [14,35]. The parallax shift problem for geostationary satellites could be treated as a special case for low-orbit satellites. This problem for low-orbit satellites could be modeled with similar equations to those presented above, however taking into account the position and orientation of the satellite in the Cartesian coordinates space.
In this paper, the parallax effect was described using an ellipsoidal earth model. However, ellipsoidal model clearly does not fully reflect the real shape of Earth. Therefore, in situations where ellipsoidal model precision is insufficient, the numerical model of Earth's gravitational field and geoid values should be utilized. In this case, it would be necessary to describe the parallax effects using differential equations and solve them using a numerical approach. In this case, the most problematic issue would be to determine perpendicular paths to the equipotential boundaries of Earth's gravitational field.