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

Abstract

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).

1. 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).

2. Nature of Parallax Shift Problem and Vicente et al./Koenig Method

2.1. 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]:

$$\{\begin{array}{c}x=\left(N\left({\phi }_g\right)+h\right)\cos {\phi }_g\cos ({\lambda }_g-{\lambda }_0)\\ y=\left(N\left({\phi }_g\right)+h\right)\cos {\phi }_g\sin ({\lambda }_g-{\lambda }_0)\\ z=\left(N\left({\phi }_g\right)\left(1-e^2\right)+h\right)\sin {\phi }_g\end{array}$$

(1)

where $N\left({\phi }_e\right)=\frac{a}{\sqrt{1-e^2{\sin }^2{\phi }_g}}$ is the prime vertical radius of the curvature, $e^2=\frac{a^2-b^2}{a^2}$ is the square of eccentricity, $a$ is Earth’s semi-major axis, $b$ is Earth’s semi-minor axis, $h$ is the cloud top height, ${\phi }_g$, ${\lambda }_g$ is the geodetic latitude and longitude, and ${\lambda }_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 ${\phi }_g$, ${\lambda }_g$ (see Figure 2). This is a more precise model than the flat-earth model or the spherical model. Note that: all longitudes (${\lambda }_g$, ${\lambda }_c$ and ${\lambda }_p$) are equal and the same. Subscripts are given to formally distinguish these values between corresponding latitudes that have different definitions (see Figure 2).

Pixel displacement in satellite view coordinates is defined as:

$$p_{disp}\left(h\right)=\sqrt{c_y^2{\left({\phi }_s\left(h\right)-{\phi }_s\left(0\right)\right)}^2+c_x^2{\left({\lambda }_s\left(h\right)-{\lambda }_s\left(0\right)\right)}^2}$$

(2)

where $c_x$ and $c_y$ are constants that allow for sensor inclination angles to be converted to pixels or distance units in the satellite view space. Also, ${\phi }_s\left(h\right)$ and ${\lambda }_s\left(h\right)$ are defined as:

$${\phi }_s\left(h\right)={\tan }^{-1}\frac{z\left(h\right)}{\sqrt{{\left(x\left(h\right)-l\right)}^2+y{\left(h\right)}^2}}$$

(3)

$${\lambda }_s\left(h\right)=-{\tan }^{-1}\frac{y\left(h\right)}{x\left(h\right)-l}$$

(4)

where $x\left(h\right)$, $y\left(h\right)$, and $z\left(h\right)$ are cloud top coordinates from Equation (1) as functions of h; $l=a+h_s$– distance from center of Earth to satellite; $a$ is the Earth’s semi-major axis; $h_s$ is the distance from the surface of Earth to the satellite. In order to illustrate $p_{disp}\left(h\right)$, 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 $c_x=c_y=\frac{h_s}{3 km/px}$). 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 ground-based radar data [7]. An example is depicted in Figure 4.

2.2. 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:

$$\{\begin{array}{c}x=R_{loc}\left({\phi }_c\right)\cos {\phi }_c\cos ({\lambda }_c-{\lambda }_0)\\ y=R_{loc}\left({\phi }_c\right)\cos {\phi }_c\sin ({\lambda }_c-{\lambda }_0)\\ z=R_{loc}\left({\phi }_c\right)\sin {\phi }_c\end{array}$$

(5)

where $a$ is Earth’s semi-major axis; $b$ is Earth’s semi-minor axis; $h$ is the cloud top height; ${\phi }_c$ and ${\lambda }_c$ are the geocentric latitude and longitude (see Figure 2), respectively; ${\lambda }_0$ is the latitude of the geostationary satellite position; and $R_{loc}\left({\phi }_c\right)$ is the local radius of ellipsoid for the geocentric latitude model:

$$R_{loc}\left({\phi }_c\right)= \frac{a}{\sqrt{{\cos }^2{\phi }_c+R_{ratio}^2{\sin }^2{\phi }_c}}$$

(6)

where:

$$R_{ratio} = \frac{a}{b}$$

(7)

Satellite position (S) is defined as:

$$\{\begin{array}{c}{x}_s=a+h_s\\ y_s=0\\ z_s=0\end{array}$$

(8)

where $a$ is Earth’s semi-major axis and $h_s$ is the distance from the surface of Earth to the satellite.

The correction procedure is as follows:

• Designate satellite position $S$ in the Cartesian coordinates system;
• Designate the position of cloud top image $I$ in the Cartesian coordinates system using Equation (5);
• Designate vector $\overrightarrow{\left|IS\right|}$;
• Designate coefficient $c$, which allows Cartesian coordinates of the cloud top to be calculated using the following equations (see Figure 5):

  $$\overrightarrow{\left|OT\right|}=\overrightarrow{\left|OI\right|}+c\overrightarrow{\left|IS\right|}$$

  (9)

  where $\overrightarrow{\left|OT\right|}$ is described by the ellipsoid parametric equation:

  $$\{\begin{array}{c}{x}_{\overrightarrow{\left|OT\right|}}=\left(a+h\right)\cos {\phi }_p\cos ({\lambda }_p-{\lambda }_0)\\ y_{\overrightarrow{\left|OT\right|}}=\left(a+h\right)\cos {\phi }_p\sin ({\lambda }_p-{\lambda }_0)\\ z_{\overrightarrow{\left|OT\right|}}=\left(b+h\right)\sin {\phi }_p\end{array}$$

  (10)

  where ${\phi }_p$ and ${\lambda }_p$ are the parametric latitude and longitude (see Figure 2). Therefore, Equation (9) can be presented as a set of equations:

  $$\{\begin{array}{c}\left(a+h\right)\cos {\phi }_p\cos ({\lambda }_p-{\lambda }_0)= x_{\overrightarrow{\left|OI\right|}}+c{x}_{\overrightarrow{\left|IS\right|}}\\ \left(a+h\right)\cos {\phi }_p\sin ({\lambda }_p-{\lambda }_0)= y_{\overrightarrow{\left|OI\right|}}+ c{y}_{\overrightarrow{\left|IS\right|}}\\ \left(b+h\right)\sin {\phi }_p= z_{\overrightarrow{\left|OI\right|}}+ c{z}_{\overrightarrow{\left|IS\right|}}\end{array}$$

  (11)

  Squaring each equation and adding them according to their sides leads to a square equation, which can be solved with respect to c:

  $$\frac{{(x_I+c{x}_{\overrightarrow{\left|IS\right|}})}^2+{(y_I+c{y}_{\overrightarrow{\left|IS\right|}})}^2}{{\left(a+h\right)}^2}+\frac{{(z_I+c{z}_{\overrightarrow{\left|IS\right|}})}^2}{{\left(b+h\right)}^2}-1=0$$

  (12)
• Apply c to calculate the Cartesian coordinates of $T$ - $x_{\overrightarrow{\left|OT\right|}}$, $y_{\overrightarrow{\left|OT\right|}}$, and $z_{\overrightarrow{\left|OT\right|}}$.
• Calculate the geocentric ellipsoidal coordinates of $T$:

  $$\{\begin{array}{c}{\phi }_c={\tan }^{-1}\frac{ z_{\overrightarrow{\left|OT\right|}}}{\sqrt{x_{\overrightarrow{\left|OT\right|}}^2+y_{\overrightarrow{\left|OT\right|}}^2}}\\ {\lambda }_c={\tan }^{-1}\frac{y_{\overrightarrow{\left|OT\right|}}}{x_{\overrightarrow{\left|OT\right|}}}+{\lambda }_0\end{array}$$

  (13)
• If required for further computation, a geodetic latitude can be calculated:

  $${\phi }_g={\tan }^{-1}\frac{a^2}{b^2}\frac{ z_{\overrightarrow{\left|OT\right|}}}{\sqrt{x_{\overrightarrow{\left|OT\right|}}^2+y_{\overrightarrow{\left|OT\right|}}^2}}$$

  (14)

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 $h$ above the ellipsoid along the normal vector. Instead, it describes the point on the ellipsoid with the semi-axes increased by $h$, therefore this method is burdened with error because of the inadequacy of the model.

3. Parallax Error Correction Methods with Lower Error

3.1. 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:

$${\phi }_g={\tan }^{-1}\frac{{\left(a+h\right)}^2}{{\left(b+h\right)}^2}\frac{ z_{\overrightarrow{\left|OT\right|}}}{\sqrt{x_{\overrightarrow{\left|OT\right|}}^2+y_{\overrightarrow{\left|OT\right|}}^2}}$$

(15)

If further computation requires the geocentric latitude, this can be calculated using the following equation:

$${\phi }_c={\tan }^{-1}\frac{b^2}{a^2}\tan {\phi }_g$$

(16)

This modification allows the correction error to be reduced to centimeters. Details will be presented in the experimental section.

3.2. 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:

$$\{\begin{array}{c}x=-q\cos {\phi }_s\cos {\lambda }_s+l\\ y=-q\cos {\phi }_s\sin {\lambda }_s\\ z=q\sin {\phi }_s\end{array}$$

(17)

where $l=a+h_s$ is the distance from Earth’s center to the satellite; $a$ is Earth’s semi-major axis; $h_s$ is the distance from the surface of Earth to the satellite; ${\phi }_s$ and ${\lambda }_s$ are satellite inclination angles; $q$ 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) should be merged, obtaining the following set of equations:

$$\{\begin{array}{c}\left(N\left({\phi }_g\right)+h\right)\cos {\phi }_g\cos ({\lambda }_g-{\lambda }_0)=-q\cos {\phi }_s\cos {\lambda }_s+l\\ \left(N\left({\phi }_g\right)+h\right)\cos {\phi }_g\sin ({\lambda }_g-{\lambda }_0)=-q\cos {\phi }_s\sin {\lambda }_s\\ \left(N\left({\phi }_g\right)\left(1-e^2\right)+h\right)\sin {\phi }_g=q\sin {\phi }_s\end{array}$$

(18)

The root of the above system of equations (${\phi }_g$ and ${\lambda }_g$) is the geodetic coordinates of point B. However, due to the entanglement of the ${\phi }_g$ 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:

$$\mathit{f}\left({\phi }_g,{\lambda }_g,q\right)=\left[f_1\left({\phi }_g,{\lambda }_g,q\right),f_2\left({\phi }_g,{\lambda }_g,q\right),f_3\left({\phi }_g,{\lambda }_g,q\right)\right]$$

(19)

where:

$$\begin{array}{c}{f}_1\left({\phi }_g,{\lambda }_g,q\right)=\left(N\right({\phi }_g)+h)\cos {\phi }_g\cos ({\lambda }_g-{\lambda }_0)+q\cos {\phi }_s\cos {\lambda }_s-l\\ f_2\left({\phi }_g,{\lambda }_g,q\right)=\left(N\right({\phi }_g)+h)\cos {\phi }_g\sin ({\lambda }_g-{\lambda }_0)+q\cos {\phi }_s\sin {\lambda }_s\\ f_3\left({\phi }_g,{\lambda }_g,q\right)=\left(N\right({\phi }_g\left)\right(1-e^2)+h)\sin {\phi }_g-q\sin {\phi }_s\end{array}$$

(20)

In Equation (19), $\left|\right|\mathit{f}\left({\phi }_g,{\lambda }_g,q\right)||$ 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:

$${\mathit{p}}_{n+1}={\mathit{p}}_n-{\left(\nabla \mathit{f}\left({\mathit{p}}_n\right)\right)}^{-1}\mathit{f}\left({\mathit{p}}_n\right)$$

(21)

where:

$$\mathit{p} :=\left[{\phi }_g,{\lambda }_g,q\right]$$

(22)

and:

$$\nabla \mathit{f}\left({\mathit{p}}_n\right)=\left[\begin{array}{c}\nabla f_1\left({\mathit{p}}_n\right)\\ \nabla f_2\left({\mathit{p}}_n\right)\\ \nabla f_3\left({\mathit{p}}_n\right)\end{array}\right]$$

(23)

and:

$$\nabla :=\left[\begin{array}{ccc}\frac{\partial }{\partial {\phi }_g}& \frac{\partial }{\partial {\lambda }_g}& \frac{\partial }{\partial q}\end{array}\right]$$

(24)

The stopping condition is defined as:

$$\left|\right|\mathit{f}\left({\mathit{p}}_n\right)\left|\right|<\epsilon$$

(25)

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):

$$\mathit{g}\left({\mathit{p}}_n\right)=\left[\begin{array}{c}{\left(f_1\left({\mathit{p}}_n\right)\right)}^2\\ {\left(f_2\left({\mathit{p}}_n\right)\right)}^2\\ {\left(f_3\left({\mathit{p}}_n\right)\right)}^2\end{array}\right]$$

(26)

with the gradient defined as:

$$\nabla \mathit{g}\left({\mathit{p}}_n\right)=\left[\begin{array}{c}2f_1\left({\mathit{p}}_n\right)\nabla f_1\left({\mathit{p}}_n\right)\\ 2f_2\left({\mathit{p}}_n\right)\nabla f_2\left({\mathit{p}}_n\right)\\ 2f_3\left({\mathit{p}}_n\right)\nabla f_3\left({\mathit{p}}_n\right)\end{array}\right]$$

(27)

and the stop condition:

$$\left|\right|\mathit{g}\left({\mathit{p}}_n\right)\left|\right|={\left|\right|\mathit{f}\left({\mathit{p}}_n\right)\left|\right|}^2<{\epsilon }^2$$

(28)

Another issue that occurs in the numerical calculation problem is the big difference in scale between ${\phi }_g, {\lambda }_g$, which are expressed in radians, and $q$, which is expressed in meters. To handle this problem, all distances ($a,b,h,l,q$) should be divided by $a$. This operation will bring $q$ to a similar scale as ${\phi }_g,{\lambda }_g$.

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.

4. 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 ${\phi }_s$ and ${\lambda }_s$ for simulated cloud top heights. With the ${\phi }_s$ and ${\lambda }_s$ 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:

• Prepare a grid of geodetic coordinates: ${\phi }_g\in \langle -90°;90°\rangle$, ${\lambda }_g\in \langle -90°;90°\rangle$, with $1°$ steps for each dimension;
• Transform the grid coordinates to the geostationary view coordinates system, ${\phi }_s$, ${\lambda }_s$ (from now on called the base ${\phi }_s$, ${\lambda }_s$) [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 (NaN – floating point special value).
• For each $h\in \left\{2 \mathrm{km}, 4 \mathrm{km}, 8 \mathrm{km}, 12 \mathrm{km}, 16 \mathrm{km}\right\}$, the following steps are performed:
  • For each ${\phi }_g$ and ${\lambda }_g$ and with $h$, calculate the $x, y, z$ coordinates using Equation (1);
  • Using $x, y, z$, calculate the geostationary view coordinates ${\phi }_s$ and ${\lambda }_s$;
  • With ${\phi }_s$, ${\lambda }_s$, and $h$, run the correction algorithms: Vicente et al./Koenig, Vicente et al./Koenig augmented, and the numerical geodetic coordinates method;
  • Each algorithm returns ${\phi }_g^{\prime }$, ${\lambda }_g^{\prime }$, which should be transformed to ${\phi }_s^{\prime }$, ${\lambda }_s^{\prime }$;
  • The distance between the simulated original base ${\phi }_s$, ${\lambda }_s$ and ${\phi }_s^{\prime }$, ${\lambda }_s^{\prime }$ 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 $h_s$ allows the result to be calculated in metric units (meters) as distances on a sphere of radius $h_s$ 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: $200m/a$;
• 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 Figure 8. The results using the geodetic coordinates numerical method are presented in Figure 9 and Figure 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.

5. 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:

$$\mathrm{dBZ}= 10lo{g}_{10}Z$$

(29)

where Z is radar reflectance. Reflectance is described as the following empirical relation with precipitation rate R [mm/h] [29]:

$$\mathrm{Z}= 200R^{1.6}$$

(30)

and cloud optical thickness (COT) in logarithmic form [30,31].

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 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.

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:

$$\frac{\partial p_{disp}\left(h\right)}{\partial h}$$

(31)

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 ${\phi }_g$ and ${\lambda }_g$. Therefore, it can be approximated as the mean slope $p_{disp}\left(h\right)$ in respect to h, for instance:

$$\frac{p_{disp}\left(12 \mathrm{km}\right)}{12 \mathrm{km}}$$

(32)

where: $c_x=c_y=h_s$. The displacement sensitivity depends on ${\phi }_g$ and ${\lambda }_g$, therefore its value varies around the globe. Sensitivity values for cities from Figure 3 are presented in

Table 1. The displacement function sensitivity in respect to h for five chosen cities for the geostationary sensor at longitude 0°. (N – North, S – South, E – East, W - West).

City	Geodetic Coordinates	Displacement Sensitivity as in Equation (32)
Cape Town	33.9253° S 18.4239° E	0.667
Madrid	40.4177° N 3.6947° W	0.696
Brasília	15.7839° S 47.9142° W	0.784
Gdańsk	54.3475° N 18.6453° E	0.827
Tromsø	69.6667° N 18.9333° E	0.868

.

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.

6. 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 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.