Downscaling of Satellite OPEMW Surface Rain Intensity Data

Abstract

This paper presents a geostatistical downscaling procedure to improve the spatial resolution of precipitation data. The kriging method with external drift has been applied to surface rain intensity (SRI) data obtained through the Operative Precipitation Estimation at Microwave Frequencies (OPEMW), which is an algorithm for rain rate retrieval based on Advanced Microwave Sounding Units (AMSU) and Microwave Humidity Sounder (MHS) observations. SRI data have been downscaled from coarse initial resolution of AMSU-B/MHS radiometers to the fine resolution of Spinning Enhanced Visible and InfraRed Imager (SEVIRI) flying on board the Meteosat Second Generation (MSG) satellite. Orographic variables, such as slope, aspect and elevation, are used as auxiliary data in kriging with external drift, together with observations from Meteosat Second Generation-Spinning Enhanced Visible and InfraRed Imager (MSG-SEVIRI) in the water vapor band (6.2 µm and 7.3 µm) and in thermal-infrared (10.8 µm and 8.7 µm). The validation is performed against measurements from a network of ground-based rain gauges in Southern Italy. It is shown that the approach provides higher accuracy with respect to ordinary kriging, given a choice of auxiliary variables that depends on precipitation type, here classified as convective or stratiform. Mean values of correlation (0.52), bias (0.91 mm/h) and root mean square error (2.38 mm/h) demonstrate an improvement by +13%, −37%, and −8%, respectively, for estimates derived by kriging with external drift with respect to the ordinary kriging.

1. Introduction

Rainfall is of primary importance in many scientific fields, such as meteorology, hydrology, agriculture, ecology and other environmental sciences [1]. Precipitation intensity can be estimated with different techniques, including rain gauge, ground-based radar, and satellite remote sensing observations. Estimates from satellite are particularly relevant, as they assure a global coverage of the Earth [2]. The major drawback of satellite rainfall remote sensing lies in their coarse spatial resolution, which hinders the investigation of spatial variability. This calls for the development of techniques to improve the spatial resolution of satellite rainfall remote sensing data.

A large number of downscaling methods has been developed and applied in the last few years. Most of these methods are based on the correlation between rainfall and environmental information such as latitude, longitude, altitude, slope, aspect and other orographic characteristics. In particular, Next Generation Radar (NEXRAD) daily precipitation fields were downscaled from 16 km to 4 km by considering orographic effects on precipitation distribution [3]. This method consists of three parts, namely the rain-pixel clustering, the multivariate regression and the random cascade. In [4], a statistical downscaling method—based on the relationships between precipitation and both terrain factors (e.g., slope, aspect and roughness) and meteorological conditions (e.g., humidity and temperature)—was proposed to disaggregate the Tropical Rainfall Measuring Mission (TRMM) 3B42 from 25 km to 1 km. Other studies introduced the positive relation between vegetation and precipitation through the use of the Normalized Difference Vegetation Index (NDVI). Among these, Ref. [5] explored the relation between TRMM rainfall estimates and NDVI at different spatial scales; the derived relation has then been used to develop a downscaling method based on an exponential regression model. In [6], a multiple linear regression model was developed using both NDVI and Digital Elevation Model (DEM) as independent variables. Ref. [7] presented a geostatistical downscaling procedure, namely a geographically-weighted regression kriging, based on the relationship between precipitation and other variables, such as NDVI and DEM, to downscale the TRMM 3B43 product from 25 km to 1 km. Two methods were used in [7]—namely a geographical difference analysis and a geographical ratio analysis—to calibrate the downscaled TRMM precipitation data. Ref. [8] applied an integrated downscaling method, based on environmental information such as vegetation, topography, drought and albedo derived from Moderate Resolution Imaging Spectroradiometer (MODIS) products. Two different downscaling approaches—a multiple linear regression and an artificial neural network—were compared and used to downscale TRMM precipitation data from 25 km to 1 km spatial resolution. The above studies demonstrate that downscaled precipitation data better capture the spatial variability compared to the original datasets.

In this work, we downscale precipitation data derived from the operational version of the Precipitation Estimation at Microwave Frequencies (PEMW), an algorithm for surface rain intensity (SRI) retrievals, developed at the Institute of Methodologies for Environmental Analysis of the National Research Council of Italy (CNR-IMAA) [9]. SRI data by PEMW is derived from Advanced Microwave Sounding Units (AMSU) and Microwave Humidity Sounder (MHS) radiometers on board Low-Earth-Orbit satellites and their spatial resolution range from 16 km at nadir to 51 km at maximum scanning angle. The Operative version of PEMW (OPEMW) was validated against simultaneous ground-based observations from weather radar systems and rain gauges [10]. OPEMW SRI data are downscaled from the AMSU-B/MHS spatial resolution (16 km at the sub satellite point—SSP) to the finer MSG-SEVIRI spatial resolution (3 km at the SSP) to possibly improve their accuracy and capture higher-resolution spatial variability. Among the methods mentioned above, a kriging technique is chosen, namely kriging with external drift, as it was shown to provide estimations and associated errors with satisfactory results even in the presence of few input data. Kriging with external drift different variables and their combinations are used as auxiliary data to improve the disaggregation process. Particularly, the Brightness Temperatures (BTs) from the MSG-SEVIRI [11] are exploited together with some environmental information, such as elevation, slope and aspect. The MSG-SEVIRI observations used here are from four channels, two centered in water vapour bands (at 6.2 µm and 7.3 µm) and two in the thermal infrared (at 10.8 µm and 8.7 µm). The precipitation episodes are divided into two types (convective and stratiform), and the trend producing the best results are searched.

The paper is organized as follows: Section 2 provides a description of the data used for the implementation and the validation of the downscaling procedure; it also describes the applied methodology, focusing on the construction of the variogram and on the choice of auxiliary variables. Section 3 shows the validation comparing downscaled surface rain rate and rain gauge measurements. In Section 4, we draw our conclusions.

2. Materials and Methods

2.1. Data Set

This section describes the OPEMW. Furthermore, we analyze the data used as auxiliary variables in kriging with external drift, such as the MSG/SEVIRI channels and the orographic variables (slope, aspect and elevation). Finally, we discuss details of the validation dataset, obtained by a network of ground-based rain gauges.

2.1.1. OPEMW

OPEMW is the operational version of PEMW running at IMAA-CNR since 2010. A detailed description of PEMW software can be found in [9]; here, we only provide a brief discussion. PEMW consists of a rainfall estimation algorithm that exploits both the radiometric observations made at window channels of 89 GHz and 150 GHz and at the water vapor band of 183 GHz [9]. It exploits the window channels to detect the size of the precipitating particles and to sense low-level precipitation, while the three water vapor bands are used to discriminate precipitation at different altitudes (high convective system, middle-altitude stratiform precipitation) [10]. PEMW is based on the observations acquired from AMSU-B on board the National Oceanic and Atmospheric Administration’s (NOAA) Polar Operational Environmental Satellites (POES), from MHS on board the European Polar System (EPS) and finally from the most recent NOAA POES. AMSU-B and MHS are respectively cross-track, line scanning microwave radiometers measuring radiances in five channels in the frequencies ranging from 89 GHz to 190 GHz. In particular, AMSU-B exploits channels at central frequencies of 89 GHz, 150 GHz, 183 ± 1 GHz, 183 ± 3 GHz and 183 ± 7 GHz, while MHS at central frequencies of 89 GHz, 157 GHz, 183 ± 1 GHz, 183 ± 3 GHz and 190 GHz. AMSU-B and MHS fly at a nominal altitude of 850 km, and they scan the Earth about ±50° from subsatellite point. Each channel has an antenna beam width of 1.1°. This provides a resolution of 16 km at nadir and 90 consecutive fields of view (FOVs) per scan. One of the main advantages related to the use of AMSU-B and MHS observations is the good spatial resolution, which is often an issue for MW instruments. The OPEMW algorithm has been validated against ground-based observations from a network of 20 weather radar systems and a network of more than 3000 rain gauges distributed over the Italian territory [10]. The data set used for the validation spans over one year of surface rain intensity (July 2011–June 2012). The validation shows 98% accuracy in correctly identifying rainy and non-rainy areas. The correlation coefficient is larger than 0.8 and 0.9 with respect to rain gauge and weather radars, respectively, though a binned analysis in the 0–15 mm ${\mathrm{h}}^{-1}$ range suggests that the algorithm tends to overestimate rain rate values below 6–7 mm ${\mathrm{h}}^{-1}$ and underestimate those above 6–7 mm ${\mathrm{h}}^{-1}$.

2.1.2. MSG-SEVIRI

SEVIRI is the visible/infrared imager on board the geostationary MSG satellite. SEVIRI is characterized by high temporal (15 min) and spatial resolutions. SEVIRI is a 50 cm diameter aperture, line by line scanning radiometer, which has the capability to observe the Earth in 12 spectral channels [11]. Out of these, eleven channels, namely the Visible (VIS), Near-InfraRed (NIR) and InfraRed (IR) ones, cover the full disk and have an imaging sampling distance of 3 km at SSP. Conversely, the remaining twelfth channel, i.e., the High Resolution Visible (HRV) channel, covers half the full disk with a 1 km at SSP imaging sampling distance at subsatellite point [12,13].

The MSG-SEVIRI BTs from the water vapor (6.2 and 7.3 µm) and the thermal-IR (10.8 and 8.7 µm) channels are exploited as auxiliary variables in the implementation of the kriging with external drift method. We consider: (i) BT differences between thermal-IR at 10.8 μm and water vapor at 6.2 μm channels, as this information is useful to identify deep-convection areas [14,15,16,17]; (ii) BT differences between water vapor at 6.2 and 7.3 μm channels, to recognize areas of intense precipitation [18]; and finally (iii) BT difference between 8.7 μm and 10.8 μm (thermal-infrared) channels, useful to discriminate liquid and ice cloud that could be associated to mid-level stratiform cloud (nimbostratus) and convective clouds, respectively.

2.1.3. Rain Gauge Network

The downscaled surface rain intensity data are validated against ground observations of rain rate measured by rain gauges distributed on the area of Basilicata region, Southern Italy (Figure 1). The considered rain gauge network is currently managed and maintained by the Functional Center Decentralized (FCD) of Basilicata [19]. The FCD carries out daily hydrometeorological monitoring and forecasting activities using meteorological data from national forecasting models and data recorded by ground stations. The FCD of Basilicata currently manages 63 monitoring stations homogeneously distributed on the regional territory. The sensors installed on each station vary depending on the monitored site and, generally, consist of rain gauges, hydrometers, thermometers, anemometers, barometers and radiometers. The data flow into an archive managed by the FCD and they are processed after quality control operations. The acquisition and processing of data measured by rain gauges are performed at different temporal intervals, ranging from 5 to 60 min. We consider only 44 rain gauges for the validation of the downscaling procedure, as these provide data at the shortest time interval (15 min). The rain rate is measured by double-tipping rain gauges and the unit of measurement of the water deposited on the ground is millimeters (mm). The rain gauges used in the validation are homogeneously distributed on the territory under analysis, as shown in Figure 1. More details on considered rain gauges are shown in Appendix B.

2.1.4. Elevation, Slope, Aspect

Digital Terrain Models (DTMs) are a primary input to any modelling or quantification process involving the earth’s topography. DTMs are raster files containing elevation data for each raster cell. DTMs are widespread for calculations, manipulations and further analysis of a geographic area, and are mainly used for extrapolating elevation information and other derived data, such as slope and aspect, related to the morphological characteristics of the study area. DTM datasets are widely used in GIS applications, and several built-in tools are available to turn the DTM into a derivative map. There are several free DTM products available, which feature high accuracy and data resolution. The accuracy of the elevation is highly correlated to orographic features and DTMs play an important role in the orographic analysis. In this work, we used a DTM from Shuttle Radar Topography Mission (SRTM) with a spatial resolution of 3 arc seconds (≈90 m) in order to assess properly the quality of our algorithm. SRTM data were collected during the 11-day mission in February 2000. Since then, SRTM data were described in detail [20,21,22] and became accessible online for free [23] SRTM3 has 90 m spatial resolution at the equator and is provided in mosaiced 5 deg × 5 deg tiles for easy download and use. All tiles are produced from a seamless dataset to allow easy mosaicking. These tiles are available in both ArcInfo ASCII and GeoTiff, formatted to ease their use in a variety of image processing and GIS applications.

The slope map was obtained using the maximum gradient of the plane tangent to the surface of the ground at a certain point. One of the methods for calculating the slope in one point is the finite differences [24], for which the following calculation example is given. The east–west and north–south gradient for the center cell in the 3 × 3 floating window (kernel) is defined as follows:

$$\left[\frac{\partial z}{\partial y}\right]=\left[\frac{\left(z_{i-1,j-1}+2\times z_{i,j-1}+z_{i+1,j-1}\right)-\left(z_{i-1,j+1}+2\times z_{i,j+1}+z_{i+1,j+1}\right)}{8\times \Delta y}\right],$$

(1)

$$\left[\frac{\partial z}{\partial x}\right]=\left[\frac{\left(z_{i+1,j+1}+2\times z_{i+1,j}+z_{i+1,j-1}\right)-\left(z_{i-1,j+1}+2\times z_{i-1,j}+z_{i-1,j-1}\right)}{8\times \Delta x}\right],$$

(2)

where z represents the elevation of the eight surrounding cells and $\Delta x,\Delta y$ specify the cell dimensions in the horizontal and lateral directions, respectively (Figure 2).

The slope angle is calculated with respect to the central cell i, j of the kernel window and is given by the surface gradient module:

$$\mathrm{slope}=\frac{180}{\mathsf{\pi }}\times \arctan \left(\sqrt{{\left[\frac{\partial z}{\partial x}\right]}^2+{\left[\frac{\partial z}{\partial y}\right]}^2}\right).$$

(3)

The slope was calculated from the SRTM3 data using the r.slope [25] algorithm within the GRASS GIS software (version 7.4.0, Free and Open Source Software). The result is a new raster data in which each pixel is associated with the slope angle, expressed in degrees (e.g., pixels with angle values equal to 0 degrees represent a flat surface).

Another product derived from the DTM is the aspect map describing the orientation—with respect to the North—of the direction of maximum slope of the plane tangent to the ground surface. This therefore represents the geographical exposure of the surface slope and it is obtained as the ratio, in a grid mesh, between the two main gradients of the grid along the x- and y-axis. As for the slope calculation, a 3 × 3 mobile window is applied to each cell (Figure 2) of the input raster and for each cell in the center of the window the aspect is calculated using an algorithm taking into account the values of the eight adjacent cells [26].

Consequently, the corresponding aspect is calculated with the following formula:

$$\mathrm{aspect}=\frac{180}{\mathsf{\pi }}\times \arctan \left(\left[\frac{\partial z}{\partial y}\right]/\left[\frac{\partial z}{\partial x}\right]\right).$$

(4)

The aspect is a spatial information indicating the direction or azimuth of a surface, and is measured in degrees in a range of values ranging from 0° to 360°. Values with 0°, 90°, 180°, and 270° indicate respectively East, North, West and South directions. The algorithm used for the aspect calculation is r.aspect of GRASS GIS software [25], which is based on the finite difference Horn method [24] using a 3 × 3 kernel.

2.2. Methodology

This section briefly presents the downscaling technique used to improve the spatial resolution of SRI data from AMSU-B/MHS spatial resolution to MSG-SEVIRI spatial resolution, i.e., the method of kriging with external drift, focusing on the construction of the empirical semivariogram and its fit. The trend used as auxiliary data to improve the downscaling process has been chosen according to the type of precipitation, classified into convective and stratiform.

2.2.1. Downscaling Technique

The ordinary kriging method is a downscaling technique to estimate a variable of interest, based on the spatial autocorrelation of data [27]. The value of a spatial process in a not observed site, i.e., $\widehat{Z}\left(s_0\right)$ with $s_0\notin \mathbf{s}=\left\{s_1{,\dots ,s}_n\right\}$, is defined as follows:

$$\widehat{Z}\left(s_0\right)=\mathbf{\lambda }\mathrm{Z}\left(\mathbf{s}\right),$$

(5)

where $\mathbf{\lambda }$ is a vector of real weights. The unbiasedness of $\widehat{Z}\left(s_0\right)$ is satisfied by its own definition and set the sum to weights to one, i.e.,

$${\displaystyle \sum }_{i=1}^n{\lambda }_i=1.$$

(6)

The method of kriging with external drift, instead, also makes use of auxiliary information. In particular, the value of $\widehat{Z}\left(s_0\right)$ is defined as follows:

$$\widehat{Z}\left(s_0\right)=\mu \left(\mathbf{s}\right)+\mathbf{\lambda }Y\left(\mathbf{s}\right),$$

(7)

where $\mu \left(\mathbf{s}\right)$ is called drift or trend and it is a linear combination of deterministic functions, i.e., $\mu \left(\mathbf{s}\right)=\mathbf{a}f\left(\mathbf{s}\right)$ with $\mathbf{a}$ vector of coefficients, $Y\left(\mathbf{s}\right)=\left(Z\left(\mathbf{s}\right)-\mu \left(\mathbf{s}\right)\right)$ represents the residual of the random process (assumed to be intrinsically stationary with mean equal to zero) and $\mathbf{\lambda }$ is a vector of the weights of the corresponding residual $Y\left(\mathbf{s}\right)$ [27].

In both cases, the weights must be determined to minimize the estimation of the variance and to ensure the estimator has no bias. The variance estimation of the differences between the value of the variable under study at two different sites $s_i$ and $s_j$ [27,28] is defined through a function called semivariogram, i.e.,:

$$\gamma \left(s_i,s_j\right)=\frac{Var\left(Z\left(s_i\right)-Z\left(s_j\right)\right)}{2}=\frac{E\left({\left|Z\left(s_i\right)-Z\left(s_j\right)\right|}^2\right)}{2}.$$

(8)

The semivariogram is a function of data variability between pairs of points at various distances and describes the probabilistic properties of a process; it is very useful in the analysis of the spatial dependence of a random process. The semivariogram is a continuous function in the origin, but, in practice, it is often observed that $\gamma \left(0\right)\ne 0$ This situation is known as nugget effect ($\gamma \left(0\right)={\tau }^2$ is called nugget) and it is linked to measurement errors or spatial resolution problems in sampling, which in turn affect the variogram value at very small distances. Another important feature of the semivariogram function is the sill, defined as:

$${\tau }^2+{\sigma }^2=\underset{\Vert h\to \infty \Vert }{\lim }\gamma \left(h\right),$$

(9)

where ${\sigma }^2$ is known as the partial sill (the sill minus the nugget) and $h$ is the distance between two different observations. In case the sill takes on a finite value, it means that the stochastic process is weakly stationary; in addition, if this occurs for a finite value of $h=h^{*}$, then $h^{*}$ is said range of the variogram. The range quantifies the distance over which two different observations can be considered correlated. After the analysis of the empirical semivariogram, a fitting model must be considered to extrapolate the spatial behaviour of the observed points to the area of interest.

In the literature, there are several theoretical semivariogram models ${\widehat{\gamma }}_i$ (e.g., linear, exponential, Gaussian, wave and circular), with known analytical properties and physical meaning of parameters. Once the empirical semivariogram function is defined and the different theoretical semivariogram models are explored, it is necessary to choose a fit criterion (Appendix A), besides the graphic inspection. These steps allow for providing the weights for the spatial interpolation through kriging techniques.

2.2.2. Choice of the Auxiliary Variables for Kriging with External Drift

A crucial step in the method of kriging with external drift lies in the selection of the auxiliary variables, which strongly depends on the parameter to be downscaled. In this work, we exploit the known relationship between rainfall and orography [1], in order to identify the most important auxiliary variables in the construction of the semivariogram. In particular, we characterize the above relationship by testing the slope, aspect and elevation data as auxiliary variables. This choice is mainly justified by the characteristics of the area under study, featuring a complex orography—i.e., a continuous alternation of mountainous areas, valleys and plains—which, in turn, affects the local climate. As precipitation may generally be classified in convective type (characterized by strong vertical velocity field, high rainfall intensity and small coverage area) and stratiform type (characterized by weak vertical velocity field, lower rainfall intensity and a more homogeneous coverage area) [29,30], we also investigated the correlation between the precipitation type and the trend to use in the method of kriging with external drift. To this aim, we took into consideration the information provided by MSG-SEVIRI, i.e., BT at four channels, in order to distinguish the cloud cover type within the AMSU/B FOV. In particular, we used:

(i)

BT differences between the 10.8 μm (thermal-infrared) and the 6.2 μm (water vapour) channels, to identify deep-convection areas [11,14,15,16,17];

(ii)

BT differences between the 7.3 μm (water vapour) and the 6.2 μm channels, to recognize areas of intense precipitation [17,18];

(iii)

BT difference between 8.7 μm and 10.8 μm (thermal-infrared) channels, useful to discriminate liquid and ice cloud that could be associated with mid-level stratiform cloud (nimbostratus) and convective clouds, respectively.

Importantly, we have firstly performed tests using one covariate at a time, and subsequently we have used combinations of different covariates as external drift to find the combination yielding the best estimate of surface rain intensity data. The validation results for the downscaled data indicate that two different trends should be used for the two analyzed rain types. In particular, the trend combinations for the convective (TREND_CR) and stratiform (TREND_SR) rain types are:

• TREND_CR = difference between 8.7 μm and 10.8 μm channels + slope + difference between 10.8 μm and 7.3 μm channels;
• TREND_SR = difference between 8.7 μm and 10.8 μm channels + slope + difference between 7.3 μm and 6.2 μm channels.

To notice that TREND_CR and TREND_SR are the function $\mu \left(\mathbf{s}\right)$ described in Formula (7). The main difference between the two trend combinations lies in the last covariate used to single out the two types of rain. The simultaneous occurrence of convective and stratiform rain in the area under study is analysed by applying both trends on the basis of the cloud mask outcome (Classification Mask Coupling of Statistical and Physics Methods, C-MACSP [14,30]), which acts as the initial selecting criterion. In detail, when C-MACSP calls for convective type rain, the trend to be used is TREND_CR; otherwise, the trend to be used is TREND_SR. We should emphasize that the combined use of orographic features (slope, aspect and elevation) and MSG-SEVIRI data, together with the separation depending on rain type, represents a novelty element within the realm of downscaling applications by means of the kriging method with external drift. In addition, the results further support this approach, since we actually find that different covariate combinations should be used to best match the rain gauge observations.

2.2.3. Validation

The downscaled SRI data were validated against rain gauge measurements. Furthermore, to show the improvement provided by the downscaling procedures compared to the initial data, the OPEMW SRI were also compared against rain gauge observations. The correlation coefficient ($corr$), the root mean square error (RMSE), the mean bias error (MBE) and the mean absolute error (MAE) are used for the validation of the results obtained by kriging with external drift, considering its different drifts. These metrics are defined as follows:

$$corr=\frac{cov\left(p,s\right)}{{\delta }_p{\delta }_s},$$

(10)

$$RMSE={\left[\frac{1}{n }{\displaystyle \sum }_{i=1}^n{(s_i-p_i)}^2\right]}^{1/2},$$

(11)

$$MBE=\frac{1}{n}{{\displaystyle \sum }}_{i=1}^n(s_i-p_i),$$

(12)

and

$$MAE=\frac{1}{n}{\displaystyle \sum }_{i=1}^n|s_i-p_i|,$$

(13)

where p denotes the rain gauges observations, s denotes the downscaled products, n is the number of observations, $cov$ indicates the covariance operator, ${\delta }_p$ is the standard deviation of p and ${\delta }_s$ is the standard deviation of s.

The validation required the space-time colocation of data from the different sources. In particular, each rain gauge rain rate has been compared with the spatially closest OPEMW value and, after the downscaling procedure, with the closest downscaled value. With regard to the temporal colocation, each OPEMW SRI value is associated with the time of the satellite overpass, which is an instantaneous observation. In order to compare the rain rate measured by OPEMW SRI with rain gauge observations, we used the 15-minute cumulative rain value, as this is the shortest time interval available for all the considered stations, and convert it as to match the rain rate units (mm h⁻¹).

3. Results

The results of the validation are presented for the ten case studies listed in

Table 1. List of case studies.

Case	Data	Satellite Overpass (UTC)	Case	Data	Satellite Overpass (UTC)
1	30 Nov 2017	06:33	6	10 Feb 2018	07:54
2	28 Dec 2017	08:24	7	13 Feb 2018	06:58
3	02 Jan 2018	08:49	8	14 Feb 2018	12:13
4	10 Jan 2018	06:54	9	20 Feb 2018	18:32
5	12 Jan 2018	06:54	10	20 Mar 2018	17:43

These cases correspond overall to 440 comparisons between observed and both original and downscaled rain data.

The continuous statistical assessment, by means of $corr$, RMSE, MBE and MAE values, is reported in

Table 2. Continuous statistics against rain gauge observations for original OPEMW data, and data downscaled with ordinary kriging (OK) and kriging with external drift (KED) methods.

Case	MAE [mm/h]	MBE [mm/h]	RMSE [mm/h]	$\mathit{c}\mathit{o}\mathit{r}\mathit{r}$
	(OPEMW, OK, KED)	(OPEMW, OK, KED)	(OPEMW, OK, KED)	(OPEMW, OK, KED)
1	(1.46, 1.39, 1.33)	(−0.72, −1.03, −0.89)	(3.08, 2.97, 2.82)	(0.62, 0.73, 0.75)
2	(0.32, 0.35, 0.34)	(−0.22, −0.20, −0.18)	(0.69, 0.67, 0.64)	(0.24, 0.32, 0.42)
3	(1.38, 1.25, 1.03)	(0.45, 0.47, 0.41)	(2.12, 1.82, 1.46)	(0.50, 0.55, 0.72)
4	(1.17, 0.70, 0.65)	(0.42, 0.60, 0.53)	(2.73, 1.22, 1.15)	(0.62, 0.72, 0.73)
5	(2.15, 1.55, 1.43)	(0.59, 0.44, 0.38)	(3.56, 2.35, 2.48)	(0.48, 0.63, 0.68)
6	(3.29, 3.09, 2.35)	(3.29, 3.09, 2.35)	(3.92, 3.63, 2.76)	(0.33, 0.33, 0.49)
7	(1.32, 0.98, 0.92)	(0.29, 0.07, 0.04)	(1.61, 1.13, 1.06)	(0.33, 0.39, 0.41)
8	(5.21, 4.99, 4.86)	(5.21, 4.99, 4.86)	(6.12, 5.75, 5.53)	(0.27, 0.31, 0.33)
9	(1.52, 1.52, 1.52)	(−1.21, −1.39, −1.39)	(2.09, 2.03, 2.01)	(0.30, 0.31, 0.34)
10	(3.86, 3.59, 3.52)	(3.30, 3.00, 2.94)	(4.61, 4.07, 3.98)	(0.31, 0.34, 0.36)

and

Table 3. Continuous statistics against rain gauge observations for original OPEMW data, and data downscaled with ordinary kriging (OK) and kriging with external drift (KED) methods.

Statistics	OPEMW	OK	KED
MAE (mm/h)	2.17	1.94	1.70
MBE (mm/h)	1.14	1.44	0.91
RMSE (mm/h)	3.05	2.58	2.38
$corr$	0.34	0.46	0.52

for each case and for the whole dataset, respectively. In detail, the continuous statistics show the comparison between rain gauge observations against the original OPEMW SRI data, the downscaled SRI data by ordinary kriging (OK) and the downscaled SRI data by kriging with external drift (KED), obtained using the trend producing the best results for each case.

The performances of the downscaling methods should be evaluated by considering also the original OPEMW SRI efficiency for the cases considered. The results show a reasonable improvement of the correlation coefficient in the comparison between rain gauge observations and downscaled data, compared to the original OPEMW SRI data. The agreement in terms of RMSE, MBE and MAE is better for downscaled data obtained by either of the two methods. In particular, the continuous statistics shows that kriging with external drift performs better than the ordinary kriging. In fact, the use of auxiliary information improves the estimate of the downscaled variable by +13% (corr), −37% (MBE), −8% (RMSE), and −12% (MAE). The computational cost to obtain estimates of rainfall data on the MSG-SEVIRI grid, which in our case consists of a number of pixels of 698, is about 1.2 s. Time refers to a Intel(R) Core(TM) i5-4460 CPU 3.20GHz, with 8 Gb RAM.

Among the case studies listed in Table 1, we show in detail two representative cases (5 and 7). In Figure 3, the MSG-SEVIRI images at VIS (0.6 µm) and IR (10.8 µm) channels and the correspondent cloud classification (C-MACSP) maps are reported for both cases. The C-MACSP maps show convection in some areas of Basilicata for case 5, while C-MACSP detects optically thick clouds in the same area for case 7. Therefore, no convection is present in the second case; however, visual inspection of the temporal sequence of MSG-SEVIRI images seem to indicate midlevel clouds—under high thick clouds detected by C-MACSP—which probably generate stratiform precipitation. The visual inspection is necessary to identify multi-layered clouds because C-MACSP does not characterize multi-layered clouds yet.

Furthermore, Figure 4 and Figure 5 show the comparison between OPEMW SRI data, rain gauge observations and downscaled data by kriging with external drift. The continuous statistic for the two cases considered is reported in

Table 4. As in Table 3 but for data from date 12 Jan 2018 at 06:54 UTC.

Statistics	OPEMW	OK	KED
MAE (mm/h)	2.15	1.55	1.43
RMSE (mm/h)	3.56	2.48	2.35
MBE (mm/h)	0.59	0.44	0.38
$corr$	0.48	0.63	0.68

and

Table 5. As in Table 3 but for data from 13 Feb 2018 at 06:58 UTC.

Statistics	OPEMW	OK	KED
MAE (mm/h)	1.32	0.98	0.92
RMSE (mm/h)	1.61	1.13	1.06
MBE (mm/h)	0.29	0.07	0.04
$corr$	0.33	0.39	0.41

.

Considering case 5 (Figure 4), the initial OPEMW SRI data set is composed by 33 pixels; out of these, only 13 feature surface rain intensity different from zero. For this case, the results show a good improvement of correlation in the comparison between rain gauge observations and downscaled data compared to the original OPEMW SRI data (0.68 against 0.48). In addition, the agreement for downscaled data in terms of RMSE and MAE is substantial (2.35 mm/h against 3.56 mm/h and 1.43 mm/h against 2.15 mm/h, respectively).

Considering case 7 (Figure 5), the initial OPEMW SRI data set consists of 34 pixels; however, only 21 of these report surface rain intensity different from zero. In this case, the results show an improvement of correlation in the comparison between rain gauge observations and downscaled data compared to the original OPEMW SRI data (0.41 against 0.33). In addition, the agreement in terms of RMSE and MAE is significant for downscaled data (1.06 mm/h against 1.61 mm/h and 0.92 mm/h against 1.32 mm/h, respectively).

In both cases, the improvement obtained by using the auxiliary variables is evident.

4. Conclusions

In this paper, we have applied a geostatistical downscaling procedure, i.e., the method of kriging with external drift, in order to improve the spatial resolution of satellite-based rainfall observations from the original resolution to finer resolution. The kriging with external drift method features several additional advantages compared to standard kriging techniques, since it relies on auxiliary variables-sampled frequently and regularly [31]—which may further minimize the error of the estimation. Therefore, a crucial step in this method lies in the selection of the auxiliary variables, which strongly depends on the parameter to be downscaled.

The known relationship between rainfall and orography [1] lead us firstly to test the slope, aspect and elevation data as auxiliary variables. Furthermore, as precipitation may generally be classified in convective and stratiform types, we also took into account the information provided by MSG-SEVIRI BT observations, acquired in the water vapor band (6.2 µm and 7.3 µm) and in thermal-IR (10.8 µm and 8.7 µm), in order to distinguish the cloud cover and precipitation types. We initially used one covariate at a time and then a combination of different covariates as external drift to find the most suitable combination yielding the best estimate for the data analysed. In detail, kriging with external drift was applied to all the proposed cases study by varying the auxiliary variables and, subsequently, analysing the validation results for each rain type separately. From this analysis, it resulted that different auxiliary variables should be used for the two rain types. In particular, we found that the trend for the convective (TREND_CR) rain cases is a combination of slope, the 8.7 µm IR channel and the difference between 10.8 μm and 7.3 μm channels, whereas, for the stratiform (TREND_SR) rain cases, the trend is a combination of the slope, the 8.7 µm IR channel and the difference between 7.3 μm and 6.2 μm channels.

To evaluate the performances of our procedure, we considered ten case studies corresponding overall to 440 comparisons between rain gauge observations and original/downscaled rain data. The statistical analysis is based on the calculation of the RMSE, MBE, MAE and correlation. The results show a reasonable improvement of the correlation coefficient in the comparison between rain gauge observations and downscaled data, compared to the original OPEMW SRI data. In addition, the agreement in terms of RMSE, MBE and MAE is better for downscaled data. In particular, the results show that kriging with external drift clearly outperforms the ordinary kriging, which only considers coordinates and distances between observations to be downscaled. In fact, the use of auxiliary information improves the estimate of the downscaled variable by +13% (corr), −37% (MBE), −8% (RMSE), and −12% (MAE). Therefore, the proposed methodology has produced results improving the statistics compared to the original OPEMW SRI data as well as the ordinary kriging. Thus, the combination of the orographic features (slope, aspect and elevation), together with MSG-SEVIRI data, represents a novelty element that allowed us to use different covariate combination for improving the quality and resolution of satellite rainfall observations.