*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{cases} \mathbf{x} = R\_{\text{loc}}(\boldsymbol{\varphi\_c}) \cos \boldsymbol{\varphi\_c} \cos(\lambda\_c - \lambda\_0) \\ \quad \boldsymbol{y} = R\_{\text{loc}}(\boldsymbol{\varphi\_c}) \cos \boldsymbol{\varphi\_c} \sin(\lambda\_c - \lambda\_0) \\ \quad \boldsymbol{z} = R\_{\text{loc}}(\boldsymbol{\varphi\_c}) \sin \boldsymbol{\varphi\_c} \end{cases} \tag{5}$$

where *a* is Earth's semi-major axis; *b* is Earth's semi-minor axis; *h* is the cloud top height; ϕ*<sup>c</sup>* and λ*<sup>c</sup>* are the geocentric latitude and longitude (see Figure 2), respectively; λ<sup>0</sup> is the latitude of the geostationary satellite position; and *Rloc*(ϕ*c*) is the local radius of ellipsoid for the geocentric latitude model:

$$R\_{\rm loc}(\varphi\_c) = \frac{a}{\sqrt{\cos^2 \varphi\_c + R\_{\rm ratio}^2 \sin^2 \varphi\_c}} \tag{6}$$

where:

$$R\_{ratio} = \frac{a}{b} \tag{7}$$

Satellite position (S) is defined as:

$$\begin{cases} x\_s = a + h\_s \\ y\_s = 0 \\ z\_s = 0 \end{cases} \tag{8}$$

where *a* is Earth's semi-major axis and *hs* is the distance from the surface of Earth to the satellite. The correction procedure is as follows:

1. Designate satellite position *S* in the Cartesian coordinates system;


$$|\overrightarrow{OT}| = |\overrightarrow{OI}| + \varepsilon |\overrightarrow{IS}|\tag{9}$$

where <sup>→</sup> |*OT*| is described by the ellipsoid parametric equation:

$$\begin{cases} \begin{aligned} \mathbf{x}\_{\stackrel{[OT]}{|OT|}} = (a+h)\cos\varphi\_{\mathcal{P}}\cos(\lambda\_{\mathcal{P}}-\lambda\_{0})\\ \mathbf{y}\_{\stackrel{\rightarrow}{|OT|}} = (a+h)\cos\varphi\_{\mathcal{P}}\sin(\lambda\_{\mathcal{P}}-\lambda\_{0})\\ \mathbf{z}\_{\stackrel{\rightarrow}{|OT|}} = (b+h)\sin\varphi\_{\mathcal{P}} \end{aligned} \end{cases} \tag{10}$$

where ϕ*<sup>p</sup>* and λ*<sup>p</sup>* are the parametric latitude and longitude (see Figure 2). Therefore, Equation (9) can be presented as a set of equations:

$$\begin{cases} (a+h)\cos\varphi\_p \cos(\lambda\_p - \lambda\_0) = \underset{|\text{OI}|}{\text{s}} + c\mathbf{x} \underset{|\text{IS}|}{\text{}} \\ (a+h)\cos\varphi\_p \sin(\lambda\_p - \lambda\_0) = \underset{|\text{OI}|}{\text{s}} + c\mathbf{y} \underset{|\text{IS}|}{\text{}} \\ (b+h)\sin\varphi\_p = z\_{\overrightarrow{|\text{OI}|}} + cz\_{\overrightarrow{|\text{IS}|}} \end{cases} \tag{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 + cx\_{\stackrel{\scriptstyle\vert IS}{\vert IS}})^2 + (y\_I + cy\_{\stackrel{\scriptstyle\vert IS}{\vert IS}})^2}{\left(a + h\right)^2} + \frac{\left(z\_I + cz\_{\stackrel{\scriptstyle\vert IS}{\vert IS}}\right)^2}{\left(b + h\right)^2} - 1 = 0\tag{12}$$


$$\begin{cases} \begin{aligned} \varphi\_{\mathcal{C}} &= \mathsf{tan}^{-1} \frac{z\_{\stackrel{\scriptstyle \mathcal{T}}} \stackrel{\scriptstyle \mathcal{T}}{\mid \partial \mathcal{T} \vert}}{\sqrt{\mathsf{x}\_{\stackrel{\scriptstyle \mathcal{T}} \mid \mathcal{T} \vert}^{2} \vert\_{\partial \mathcal{T} \vert}}} \end{aligned} \\ \begin{aligned} \begin{aligned} \lambda\_{\mathcal{C}} &= \mathsf{tan}^{-1} \frac{\mathsf{Y}\_{\stackrel{\scriptstyle \mathcal{T}}} \stackrel{\scriptstyle \mathcal{T}}{\mid \partial \mathcal{T} \vert}}{\mathsf{x}\_{\stackrel{\scriptstyle \mathcal{T}} \mid \partial \mathcal{T} \vert}} + \lambda\_{0} \end{aligned} \end{cases} \end{cases} \tag{13}$$

7. If required for further computation, a geodetic latitude can be calculated:

$$\varphi\_{\mathcal{S}} = \tan^{-1} \frac{a^2}{b^2} \frac{z\_{|\vec{OT}|}}{\sqrt{\chi^2\_{|\vec{OT}|} + y\_{|\vec{OT}|}^2}} \tag{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.

**Figure 5.** Vector notation in the Vicente et al./Koenig method.

#### **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:

$$q\_{\vec{\mathcal{S}}} = \tan^{-1} \frac{(a+h)^2}{\left(b+h\right)^2} \frac{z\_{\vec{\mathcal{O}}|\vec{\mathcal{T}}|}}{\sqrt{x\_{\vec{\mathcal{O}}|\vec{\mathcal{T}}|}^2 + y\_{\vec{\mathcal{O}}|\vec{\mathcal{T}}|}^2}} \tag{15}$$

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

$$
\varphi\_{\varepsilon} = \tan^{-1} \frac{b^2}{a^2} \tan \varphi\_{\mathcal{S}} \tag{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{cases} \begin{cases} \begin{aligned} \mathbf{x} = -q \cos \varphi\_{\mathbf{s}} \cos \lambda\_{\mathbf{s}} + l \\ y = -q \cos \varphi\_{\mathbf{s}} \sin \lambda\_{\mathbf{s}} \\ z = q \sin \varphi\_{\mathbf{s}} \end{aligned} \end{cases} \tag{17}$$

where *l* = *a* + *hs* is the distance from Earth's center to the satellite; *a* is Earth's semi-major axis; *hs* is the distance from the surface of Earth to the satellite; ϕ*<sup>s</sup>* and λ*<sup>s</sup>* 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{cases} \begin{pmatrix} \mathcal{N}(\boldsymbol{\varrho}\_{\mathcal{S}}) + h \right) \cos \boldsymbol{\varrho}\_{\mathcal{S}} \cos(\boldsymbol{\lambda}\_{\mathcal{S}} - \boldsymbol{\lambda}\_{0}) = -q \cos \boldsymbol{\varrho}\_{s} \cos \boldsymbol{\lambda}\_{s} + l\\ \begin{pmatrix} \mathcal{N}(\boldsymbol{\varrho}\_{\mathcal{S}}) + h \end{pmatrix} \cos \boldsymbol{\varrho}\_{\mathcal{S}} \sin(\boldsymbol{\lambda}\_{\mathcal{S}} - \boldsymbol{\lambda}\_{0}) = -q \cos \boldsymbol{\varrho}\_{s} \sin \boldsymbol{\lambda}\_{s} \end{cases} \tag{18}$$
 
$$\begin{cases} \mathcal{N}(\boldsymbol{\varrho}\_{\mathcal{S}}) \left(1 - c^{2} \right) + h \right) \sin \boldsymbol{\varrho}\_{\mathcal{S}} = q \sin \boldsymbol{\varrho}\_{s} \end{cases} \tag{18}$$

The root of the above system of equations (ϕ*<sup>g</sup>* and λ*g*) is the geodetic coordinates of point *B*. However, due to the entanglement of the ϕ*<sup>g</sup>* 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:

$$f(\boldsymbol{\varphi}\_{\mathcal{S}'}, \boldsymbol{\lambda}\_{\mathcal{S}'}, q) = \left[ f\_1(\boldsymbol{\varphi}\_{\mathcal{S}'}, \boldsymbol{\lambda}\_{\mathcal{S}'}, q)\_{\prime} f\_2(\boldsymbol{\varphi}\_{\mathcal{S}'}, \boldsymbol{\lambda}\_{\mathcal{S}'}, q)\_{\prime} f\_3(\boldsymbol{\varphi}\_{\mathcal{S}'}, \boldsymbol{\lambda}\_{\mathcal{S}'}, q) \right] \tag{19}$$

where:

$$\begin{array}{l} f\_1(\boldsymbol{\varrho}\_{\mathcal{S}'}, \boldsymbol{\lambda}\_{\mathcal{S}'}, q) = \left( \mathcal{N}(\boldsymbol{\varrho}\_{\mathcal{S}}) + h \right) \cos \boldsymbol{\varrho}\_{\mathcal{S}} \cos \left( \boldsymbol{\lambda}\_{\mathcal{S}} - \boldsymbol{\lambda}\_0 \right) + q \cos \boldsymbol{\varrho}\_{\mathcal{S}} \cos \boldsymbol{\lambda}\_{\mathcal{S}} - h \\\ f\_2(\boldsymbol{\varrho}\_{\mathcal{S}'}, \boldsymbol{\lambda}\_{\mathcal{S}'}, q) = \left( \mathcal{N}(\boldsymbol{\varrho}\_{\mathcal{S}}) + h \right) \cos \boldsymbol{\varrho}\_{\mathcal{S}} \sin \left( \boldsymbol{\lambda}\_{\mathcal{S}} - \boldsymbol{\lambda}\_0 \right) + q \cos \boldsymbol{\varrho}\_{\mathcal{S}} \sin \boldsymbol{\lambda}\_{\mathcal{S}} \\\ f\_3(\boldsymbol{\varrho}\_{\mathcal{S}'}, \boldsymbol{\lambda}\_{\mathcal{S}'}, q) = \left( \mathcal{N}(\boldsymbol{\varrho}\_{\mathcal{S}}) \left( 1 - e^2 \right) + h \right) \sin \boldsymbol{\varrho}\_{\mathcal{S}} - q \sin \boldsymbol{\varrho}\_{\mathcal{S}} \end{array} \tag{20}$$

In Equation (19), *f* ϕ*g*, λ*g*, *q*  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:

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

where:

$$\mathfrak{p} := \left[ \varphi\_{\mathbb{S}'} \lambda\_{\mathbb{S}'} \mathfrak{q} \right] \tag{22}$$

and:

$$\nabla f(\boldsymbol{p}\_n) = \begin{bmatrix} \nabla f\_1(\boldsymbol{p}\_n) \\ \nabla f\_2(\boldsymbol{p}\_n) \\ \nabla f\_3(\boldsymbol{p}\_n) \end{bmatrix} \tag{23}$$

and:

$$
\nabla := \begin{bmatrix}
\frac{\partial}{\partial \psi\_{\mathcal{S}}} & \frac{\partial}{\partial \lambda\_{\mathcal{S}}} & \frac{\partial}{\partial q}
\end{bmatrix} \tag{24}
$$

The stopping condition is defined as:

$$\left\| f(p\_n) \right\| < \varepsilon \tag{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):

$$\mathbf{g}(p\_n) = \begin{bmatrix} \left(f\_1(p\_n)\right)^2\\ \left(f\_2(p\_n)\right)^2\\ \left(f\_3(p\_n)\right)^2 \end{bmatrix} \tag{26}$$

with the gradient defined as:

$$\nabla \mathbf{g} \begin{pmatrix} p\_n \\ \end{pmatrix} = \begin{bmatrix} 2f\_1(p\_n) \nabla f\_1(p\_n) \\ 2f\_2(p\_n) \nabla f\_2(p\_n) \\ 2f\_3(p\_n) \nabla f\_3(p\_n) \end{bmatrix} \tag{27}$$

and the stop condition:

$$\left\| \left[ \mathbf{g}(\boldsymbol{\nu}\_{n}) \right] \right\| = \left\| \boldsymbol{\mathfrak{f}}(\boldsymbol{\nu}\_{n}) \right\|^{2} < \varepsilon^{2} \tag{28}$$

Another issue that occurs in the numerical calculation problem is the big difference in scale between ϕ*g*, λ*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 ϕ*g*, λ*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.

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

### **4. Parallax E**ff**ect 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 ϕ*<sup>s</sup>* and λ*<sup>s</sup>* for simulated cloud top heights. With the ϕ*<sup>s</sup>* and λ*<sup>s</sup>* 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:

	- a. For each ϕ*<sup>g</sup>* and λ*<sup>g</sup>* and with *h*, calculate the *x*, *y*, *z* coordinates using Equation (1);
	- b. Using *x*, *y*, *z*, calculate the geostationary view coordinates ϕ*<sup>s</sup>* and λ*s*;
	- c. With ϕ*s*, λ*s*, and *h*, run the correction algorithms: Vicente et al./Koenig, Vicente et al./Koenig augmented, and the numerical geodetic coordinates method;
	- d. Each algorithm returns ϕ *<sup>g</sup>*, λ *<sup>g</sup>*, which should be transformed to ϕ *<sup>s</sup>*, λ *s*;

e. The distance between the simulated original base ϕ*s*, λ*<sup>s</sup>* and ϕ *<sup>s</sup>*, λ *<sup>s</sup>* 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 *hs* allows the result to be calculated in metric units (meters) as distances on a sphere of radius *hs* 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:


The results of the simulation using the Vicente et al./Koenig method and its augmented version are presented in Figures 7 and 8. The results using the geodetic coordinates numerical method are presented in Figures 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 103 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 103. 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.

**Figure 7.** Maps showing error for the Vicente et al./Koenig (**a**,**c**,**e**,**g**,**i**) method and its augmented version (**b**,**d**,**f**,**h**,**j**) for several chosen cloud heights. Error are given in meters for the geostationary satellite coordinate system. NaN values for in-scope regions occur where the algorithm failed to calculate a solution. For each map, the median (med.) error was calculated.

**Figure 8.** Error histograms for the Vicente et al./Koenig method (**a**,**c**,**e**,**g**,**i**) and its augmented version (**b**,**d**,**f**,**h**,**j**) for several chosen cloud heights. The Y-axis represents a count of 1 degree pixels, and the X-axis is the error in meters for the geostationary satellite coordinate system.

**Figure 9.** Maps with error (**a**,**c**,**e**,**g**,**i**) and number of iterations (**b**,**d**,**f**,**h**,**j**) for geodetic coordinates numerical algorithm (Geod. num. alg.) for several chosen cloud heights. Errors are given in meters for the geostationary satellite coordinate system. For each case, the median (med.) error was calculated.

**Figure 10.** Histograms of error (**a**,**c**,**e**,**g**,**i**) and number of iterations (**b**,**d**,**f**,**h**,**j**) for geodetic coordinates algorithm (Geod. num. alg.) method for several chosen cloud heights. Error is given in meters for the geostationary satellite coordinate system.

#### **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:

$$\text{dBZ} = \ 10 \log\_{10} \text{Z} \tag{29}$$

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

$$Z = \, 200 \text{R}^{1.6} \tag{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 Figures 6 and 12 still differ. In this context, it is important to note that these differences are caused by other factors that influence data acquisition, namely:


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

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\_{\rm disp}(h)}{\partial h} \tag{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 ϕ*<sup>g</sup>* and λ*g*. Therefore, it can be approximated as the mean slope *pdisp*(*h*) in respect to *h*, for instance:

$$\frac{p\_{\rm disp} \,(12 \text{ km})}{12 \text{ km}} \tag{32}$$

where: *cx* = *cy* = *hs*. The displacement sensitivity depends on ϕ*<sup>g</sup>* and λ*g*, therefore its value varies around the globe. Sensitivity values for cities from Figure 3 are presented in Table 1.

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


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.

**Author Contributions:** Conceptualization, T.B.; methodology, T.B.; software, T.B.; validation, T.B.; formal analysis, T.B.; investigation, T.B.; resources, T.B.; data curation, T.B.; writing—original draft preparation.; writing—review and editing, T.B.; visualization, T.B. All authors have read and agreed to the published version of the manuscript.

**Funding:** The research was supported under ministry subsidy for research for Gdansk University of Technology.

**Acknowledgments:** I would like to thank Andrzej Chybicki, PhD, and Tomasz Berezowski, PhD, for scientific and editorial support, as well as Marek Moszy ´nski for supervising my work. Calculations were carried out thanks to the Academic Computer Centre in Gda ´nsk. Meteorological data from on-ground radars were provided by the Polish Institute of Meteorology and Water Management, National Research Institute.

**Conflicts of Interest:** The author declares no conflict of interest.
