Geographic segmentation, spatial dependencies, and evaluation of the relative position of rain-gauges based on gridded data of mean monthly precipitation: application in Nigeria

Over the last years, high resolution grids of various climatic parameters have been developed at global scale (e.g., Hijmans et al. 2005; Mitchell & Jones 2005; Sheffield et al. 2006) providing significant information for regions with no meteorological stations. The availability of such data stimulates interest for the implementation of various statistical techniques, which have already been used to capture the spatial and seasonal patterns of precipitation or other climatic parameters using data directly from meteorological stations (Wackernagel 2003; Chen et al. 2009; Buttafuoco et al. 2011; Marchetti et al. 2015).

The development of gridded climatic data is usually performed either using climatic models (e.g., general circulation models) (Sheffield et al. 2006; Shaffrey et al. 2009; Watanabe et al. 2010) in combination with downscaling techniques (Wilby & Wigley 1997), or using interpolation techniques (Ninyerola et al. 2000; Boer et al. 2001; Hijmans et al. 2005). Advanced interpolation techniques take into account not only the site-specific observations from stations and their location but also the effects of topography and other parameters (Ninyerola et al. 2000; Hijmans et al. 2005), while climatic models include also the effects of ocean and air masses’ circulation (Sheffield et al. 2006; Shaffrey et al. 2009; Watanabe et al. 2010).

Gridded multi-temporal climatic datasets can be used in order to define regions of similar spatiotemporal variability based on one or more climatic parameters using techniques such as cluster analysis in GIS environment. Such techniques have been used on land surface temperature (LST) gridded data in an attempt to define sub-regions with different seasonal LST variability, to assess its sensitivity to climatic change and to support environmental analysis (Miliaresis 2009, 2012; Miliaresis & Partsinevelos 2010; Maeda & Hurskainen 2014). In addition, Miliaresis (2013, 2014) proposed a method for standardizing multi-temporal LST imagery in order to reveal and describe thermal anomalies using elevation, latitude and longitude dependencies using correlation and principal component analysis (PCA). These techniques have also been adopted in order to analyze the spatiotemporal patterns of precipitation and reference crop evapotranspiration allowing the terrain segmentation of a territory based on the spatiotemporal variation of the two climatic parameters (Demertzi et al. 2014; Aschonitis et al. 2016).

The climatic regionalization/segmentation of a territory based on the spatiotemporal variation of a gridded climatic parameter and its dependency on spatial features (e.g., elevation, latitude, longitude, etc.) can further be used to evaluate the representativeness of existing meteorological stations to describe the climatic variability of adjacent territories with no observed data. This information could be used as an additional tool for the evaluation of existing networks of meteorological stations, especially for rain-gauge stations, since precipitation is the most variable and one of the most important climatic parameters. The design of the spatial structure/distribution of rain-gauge networks is very complex and efforts regarding its evaluation have been performed using satellite images (Caselles & Melia 1985), radar observations (Abdella & Alfredsen 2010), geostatistical methods and optimization algorithms (Filippini et al. 1994; Papamichail & Metaxa 1996; Tsintikidis et al. 2002; Nour et al. 2006; Adhikary et al. 2015), probabilistic approaches (Shafiei et al. 2014), ‘cost–benefit’ approaches using a step-by-step reduction of stations’ density (Vose & Menne 2004; Horel & Dong 2010), entropy approaches (Krstanovic & Singh 1992; Al-Zahrani & Husain 1998; Kawachi et al. 2001; Chen et al. 2008; Yoo et al. 2008), approaches based on intensity–duration–frequency (Chebbi et al. 2013), and combinations of the aforementioned approaches (Pardo-Igúzquiza 1998; Barca et al. 2008; Cheng et al. 2008). Considering the above, terrain segmentation/clustering approaches based on available gridded data could assist the aforementioned procedures for the evaluation and design of rain-gauge networks.

The aim of the study is to present a combination of techniques for spatiotemporal analysis of multi-temporal gridded climatic datasets, which can further be used to evaluate the relative position of existing meteorological stations. The techniques were applied for the analysis of gridded mean monthly precipitation data for the period 1950–2000 and for the evaluation of the existing rain-gauge network of Nigeria.

The study area is the country of Nigeria which is confined between 4 ° to 14 ° north and 2 ° to 15 ° east covering a total area of ∼923,770 km². The selection of the specific study area was based on the fact that the climate of the country presents extremely high latitudinal variability, which ranges from south to north as follows: tropical monsoon climate (Am), tropical savanna climate (Aw), hot semi-arid climate (BSh), hot desert climate (BWh) (Peel et al. 2007) (Figure S1(a) in Supplementary material, available with the online version of this paper). Nigeria has only two seasons, a dry and a rainy/monsoon season. The dry season is influenced by dry north-easterly winds, which regulate the flow of dust-laden air masses from the Sahara Desert, known as Continental Tropical (cT) air mass (locally known as Harmattan). The monsoon season is influenced by moist south-westerly winds that regulate the flow of air masses known as the Maritime Tropical (mT) air mass which reaches the country from the Gulf of Guinea coast (tropical Atlantic). The moving boundary between the two wind systems creates a zone of moisture discontinuity known as the Inter Tropical Discontinuity (ITD) zone. mT is responsible for Nigeria's rainy season, invades the country usually after February from the south and reaches the northern part of Nigeria during June. mT invasion is followed by the northward retreat of cT and both of them control the northward shift of ITD. Precipitation variability in Nigeria is not only regulated by the ITD mechanism but also by other tropical and extra-tropical factors, such as the tropical easterly jet (TEJ), the sea surface temperature (SST) anomaly, the stratospheric Quasi-Biennial Oscillation (QBO), the El Nino Southern Oscillation (ENSO), Southern Oscillation Index (SOI) and North Atlantic Oscillation (NAO), and the Atlantic Multi-Decadal Oscillation (AMO) (Olaniran 1983, 1988a, 1988b, 1990, 1991a, 1991b; Adefolalu & Oguntoyinbo 1985; Adefolalu 1988; Olaniran & Sumner 1989a, 1989b, 1990; Bello 1996, 1998; Odekunle 2010; Oguntunde et al. 2011, 2014; Alli et al. 2012; Ogungbenro & Morakinyo 2014).

Gridded precipitation data for Nigeria were obtained from the WorldClim database (http://www.worldclim.org/). The database provides mean monthly values of precipitation at 30 arc-sec spatial resolution (∼1 × 1 km) for the period 1950–2000 at global scale (Hijmans et al. 2005). The dataset was produced using the thin-plate smoothing spline algorithm implemented in the ANUSPLIN package for interpolation, using latitude, longitude, and elevation as independent variables and mean monthly precipitation data from rain-gauge stations (Hijmans et al. 2005). The authors performed: (a) uncertainty quantification arising from the input data and the interpolation by mapping weather station density; (b) analysis of elevation bias in the weather stations and of elevation variation within grid cells; and (c) data partitioning and cross validation. A revised version of the GTOPO30 DEM based on SRTM DEM at 30 arc-sec spatial resolution was also obtained by the WorldClim database.

Location data of the current network of 33 rain-gauges in Nigeria were obtained by NOAA (www.nws.noaa.gov/tg/siteloc.php). The position of rain-gauge stations in comparison to (a) Köppen climate classification (Peel et al. 2007), (b) different ecoregions (Olson et al. 2001) and hydrologic network, and (c) elevation are given in Figure S1(a)–S1(c), respectively.

A combination of three statistical methods consisting of correlation, principal components (PCs), and cluster analysis was used in GIS environment (Arc-Info 9.3, ESRI-GIS) for imagery analysis. Correlation analysis and PCA were used in order to investigate the spatiotemporal dependencies of P according to the proposed procedures by Miliaresis (2012) while cluster analysis was used to perform terrain segmentation based on the spatiotemporal variation of mean monthly precipitation.

The correlation analysis derived the Pearson correlation coefficient R of mean monthly precipitation P versus elevation H, latitude Lat, and longitude Lon. The analysis was performed using the following toolbox in Arc-Map (Spatial Analyst>Multivariate>Band Collection Statistics) which produces a text file with the correlations’ table.

Cluster analysis was used to segment Nigeria into a number of distinct territories (clusters) with similar intensity and temporal variability of mean monthly precipitation. Cluster signatures were derived by the Iso Cluster algorithm (Ball & Hall 1965), which uses a modified iterative optimization procedure known as the migrating means technique. The algorithm separates all cells into a predefined number of distinct unimodal groups (clusters) in the multidimensional space of the input bands. The algorithm iteratively computes the minimum Euclidean distance when assigning each candidate cell to a cluster. The process starts with arbitrary means assigned by the software for each cluster (the number of clusters is given by the user). Every cell is assigned to the closest of these means (all in the multidimensional attribute space). New means are iteratively recalculated for each cluster based on the attribute distances of the cells that belong to the cluster (200 iterations with sampling interval of 10 pixels were used in this study). The optimal number of clusters to specify is usually unknown. For this reason, an initial trial is performed using a high number in order to allow the algorithm to give its maximum number of stable clusters, which is smaller than the specified one by the user (Ball & Hall 1965; Richards 1986). The user can choose the maximum number of stable clusters or to reduce their number by repeating the procedure. In our case, the maximum number was used while the value of 50 km² was assigned as the minimum threshold of total cluster area coverage. The aforementioned procedure was performed using the following toolbox in Arc-Map (Spatial Analyst>Multivariate>Iso Cluster) which produces a signatures file of variances and covariances of each cluster. The maximum likelihood classification (MLC) was used for terrain segmentation using the cluster signatures from the Iso Cluster. The MLC algorithm is based on Bayes’ theorem of decision (Gelman et al. 2014) considering that the cells in each cluster sample in the multidimensional space are normally distributed. The MLC considers both the variances and co-variances of the cluster signatures when assigning each cell to one of them. Each cluster is characterized by a mean vector and a covariance matrix. Taking into account these two characteristics, statistical probabilities are computed for the pixels to determine their proximity to the mean characteristics of the clusters to which they belong. Thus, MLC produced in our case a raster which shows the levels of % confidence of classification reliability for each classified pixel (eight reliability classes were used). A 0% reject fraction was used in order not to allow unclassified pixels. The aforementioned procedure was performed using the following toolbox in Arc-Map (Spatial Analyst>Multivariate>MLC) which produces two rasters: the raster of clusters classification and the raster of reliability classification classes.

Finally, the distribution of the 33 existing rain-gauges’ network was evaluated based on (a) the number of rain-gauges which exist in each cluster and (b) the reliability class values of the pixels in the positions of rain-gauges.

The 30 arc-sec resolution (∼1 × 1 km) maps of elevation H and mean annual precipitation P for the period 1950–2000 are given in Figure 1(a) and 1(b), respectively, while the respective cumulative area (based on the number of pixels) for the observations of elevation H and P in Figure 1(a) and 1(b) are given, respectively, in Figure 1(c) and 1(d). The respective overall mean annual value of precipitation and elevation over Nigeria are 1,174 ± 547 mm yr⁻¹ and 329 ± 222 m above the mean sea level (m.a.s.l).

The rainfall map (Figure 1(b)) shows that the highest annual P values are observed along the southern coastline while the lowest annual P values were observed in the northern part of the country. These trends indicate a clear decrease in rainfall amount with the increase of latitude verifying the results of previous studies (Bello 1998; Nicholson et al. 2000; Oguntunde et al. 2011). Taking into account the 1,077,730 pixel values of mean annual P raster, a simple regression was performed with latitude values in order to assess the latitudinal effect on annual P variability. The best fit was performed by a reciprocal-Lat function of P (P= 16,285.1/Lat ‒ 641.9, p < 0.001, R² = 0.89) (Figure S2(a) in Supplementary material, available with the online version of this paper). The regression showed that the slope, which describes the decrease of annual precipitation per latitude unit increase, is steeper in the low latitude areas. The mean latitudinal effect for the entire country was found equal to 212 mm of P reduction per 1 ° of latitude increase. Splitting Nigeria into two regions between 4 °–9 ° N and 9 °–14 ° N the respective rates were found to be 468 and 130 mm of P reduction per 1 ° of latitude increase. The graphs between mean annual P versus longitude Lon and elevation H are also given in Figure S2(b) and S2(c), respectively, without providing fitting analysis due to weak correlation coefficients.

The 30 arc-sec resolution (∼1 × 1 km) maps of mean annual precipitation P for the period 1950–2000 showed that the overall mean annual value of precipitation over Nigeria is 1,174 ± 547 mm yr⁻¹. Oguntunde et al. (2011), who also analyzed other gridded precipitation data of coarser resolution (0.5 °, CRU TS 2.1 dataset) (Mitchell & Jones 2005) for the period 1901–2000, found that the mean ± st. dev. precipitation is equal to 1,170 ± 109 mm yr⁻¹. Although the mean annual precipitations given by the two different datasets are in proximity, the difference in their standard deviations is extremely high. The significantly lower standard deviation of the second dataset suggests a significant loss of detail in local precipitation patterns due to the lower resolution. This loss of detail in the general climatology of the country has also been identified by Hassan et al. (2009).

The minimum, maximum, and mean ± st. dev. values of monthly precipitation for the studied period are also presented in Figure 2(a). According to Figure 2(a), the minimum and mean monthly precipitations show a unimodal response while the maximum monthly precipitation shows a bimodal response with two maximum extremes during June (641 mm) and September (527 mm).

The correlations of the mean monthly values of precipitation P versus elevation H, latitude Lat, and longitude Lon are given in Figure 2(b). It is indicative that the monthly variation of correlations of P versus H, Lat, and Lon follow similar patterns. All correlations are negative showing a relative stability during the period of October–June. The correlations’ magnitude of P versus these three parameters shows the following order Lat > H > Lon during this period. On the other hand, the period of July–September indicates a climatic disturbance, mainly due to ITD northward movement, which forces P-H, P-Lat, P-Lon relations to shift towards positive values. This disturbance succeeds in inverting the effects of H and Lon on monthly P variation in August (positive R) but not the effect of Lat, although it shows more intense change in comparison to H and Lon. This shift is also responsible for the bimodal precipitation patterns in the southern regions during the wet season. It is indicative that the monthly variation of absolute correlation for P-Lat is maximized during the onset (April–May) and cessation (October) of rainfall reaching values |R| > 0.9. The negative correlation between monthly P and elevation H is probably related to the fact that the higher altitude areas are located in the drier northern regions.

In order to further explore the monthly correlations observed in Figure 2(b), PCA was performed and the corresponding eigenvalues and factor loadings are presented in Table S1 (available with the online version of this paper). The first two PCs (PC-1 and PC-2) account for 91.7% of the total variance in the mean monthly P and they are given in Figure 3. The PC-1 map (Figure 3(a)) explains 80.3% of the total variance in the mean monthly P and amplifies the difference between August and the rest of the period of September–July (Table S1 and Figure 3(c)). The factor loadings of PC-1 (Figure 3(c)) are relatively stable during the period of September–July while during August they are significantly reduced. PC-1 is spatially maximized in the southern tropical zone while it is minimized in northern Nigeria, providing a first indication of which regions are more likely to present higher amounts of rainfall. PC-2 (Figure 3(b)) explains 11.4% of the total variance in the mean monthly P and amplifies the difference between the periods of October–May and June–September. The difference between these two periods is spatially maximized in the central mountainous regions of Plateau and Kaduna provinces and in the southeastern regions of the coastline. PC-2 is minimized in the southwestern region (rain forests) of Nigeria and in the extreme northeastern region (Figure 3(b)). PC-1 and PC-2 (Figure 3) provide a good visualization of the northward movement effects of IDT on precipitation patterns and the period of occurrence (period of June–September with a peak in August according to Figure 3(c)).

Cluster analysis revealed 47 clusters and their characteristics such as percent area coverage, mean elevation and mean monthly P values are given in Table 1 while their spatial distribution is given in Figure 4(a). The levels of % confidence of classification reliability for each classified pixel are also given in Figure 4(b).

The clusters were further grouped based on P modality. The clusters 1 up to 18 present unimodal while the clusters 19 up to 47 present bimodal monthly variation of P. The clusters belonging in the bimodal group present higher precipitation in comparison to the unimodal group (Table 1) justifying the bimodal response of the maximum monthly P variation for the whole country (Figure 2(a)). The spatial threshold between the unimodal and bimodal responses of P is given in Figure 5. The % coverage of the unimodal and bimodal clusters was found equal to 69% and 31%, respectively (Table 1), justifying the dominance of the first group in the unimodal response of the mean monthly P variation for the whole country (Figure 2(a)). The general characteristics of the modality groups are the following:

• The unimodal group is divided into two subgroups where the first (code 1.1) presents a unique P maximum extreme during August and the second one (code 1.2) during September (Table 1, Figure 5).

• The bimodal group is divided into three subgroups where: (a) for the first subgroup (code 2.1), the minimum extreme, which is responsible for the bimodal rainy season, appears during July while September is the most rainy month; (b) for the second subgroup (code 2.2), the minimum extreme, which is responsible for the bimodal rainy season, appears during August while September is the most rainy month; and (c) for the third subgroup (code 2.3), the minimum extreme, which is responsible for the bimodal rainy season, appears during August while June or July is the most rainy month (Table 1, Figure 5).

• The clusters 6, 9, 10, 11, 14, 16, and 42 present a distinct geographical dispersion in two different geographical regions (e.g., see black arrows for cluster 14 in Figure 4(a)). These clusters (except 42) are on the west and the east side of the upland regions of Bauchi, Kaduna, Kano, and Plateau states (located in the clusters 7, 8, 12, 13). This dispersion probably occurs because these upland regions force part of the incoming vapor masses from the south through mT to follow two different pathways through Niger and Benue River valleys (Figure S1(a)), respectively.

• The mean monthly values of P for each cluster (Table 1) also provide the length of the rainy season and the total rainfall during the rainy season over the different areas of the country. The clusters 46 and 47 present the higher annual precipitation exceeding the threshold of 3,000 mm/yr while clusters 1 and 2 present the lower annual precipitation (below the threshold of 500 mm/yr). The clusters 1–13 and 16 present at least one month with no precipitation while the highest mean monthly precipitation appears in cluster 47 during June with a value of 541 mm/month.

The derived terrain segmentation of Figure 4(a) was used to evaluate the position of rain-gauges based on the cluster characteristics. Analyzing the position of the 33 existing stations (Figure 5), it was found that 22 out of 47 clusters do not contain a rain-gauge, while 17 and 8 clusters contain one and two stations, respectively (Table 1). The indication of reliability class <5% showed that 13 stations are located in regions of transitional precipitation patterns. In the case of clusters with more than one station, the reliability class can provide a good criterion for selecting the most representative one.

A coarser clustering, as is described by the zones of Figure 5, can also be used to analyze stations’ distribution. In this case, it was observed that the number of stations that are located in the modality subgroup zones with codes 1.1, 1.2, 2.1, 2.2, and 2.3 were 12, 3, 3, 12, and 3, respectively (Table 2). Their distribution is generally homogeneous in the three zones with codes 1.1, 2.1, and 2.2. In the case of 1.2 zone, the number of stations per unit area is the lowest and the stations are located close to the center of the zone and not at its west and east side where Niger and Benue River valleys are located. For the case of zone 2.3, it was observed that the southeast side contains only one station (Figure 5). In zone 2.3 are included the clusters 46 and 47 (Figure 4(a)) which are described by the highest annual precipitation (Table 1) without having a station within their boundaries.

The application of the proposed methods using gridded climatic data for the evaluation of meteorological station networks presents the following limitations:

• (a) The minimum period of temporal resolution of available gridded climatic data is usually a month. This attribute is a limiting factor, especially in the case of precipitation, since the proposed methods cannot capture differences in frequency, duration, and intensity of rainfall. The aforementioned limitations in the case of rainfall can lead to geographical dispersion of clusters (as in the case of clusters 6, 9, 10, 11, 14, 16, and 42). Clear dispersion of clusters (e.g., cluster 14, see black arrows in Figure 4(a)) can further be used to create sub-clusters. In the case of other climatic parameters (e.g., temperature, solar radiation, etc.), the proposed methodology is expected to be more robust.

• (b) The results of the proposed methodology are highly dependent on the methods and data used for grid development. One of the most important problems in climatic grid development is the high uncertainty of grid values in high altitude regions or vast inaccessible regions where meteorological stations usually do not exist (Hijmans et al. 2005).

Due to the aforementioned limitations of the proposed methodology for evaluating the station network, it is recommended to be used as a supporting tool in combination with other methods such as those mentioned in the Introduction section. Some suggestions for future applications are:

• (a) use of the proposed methodology as a primary analysis for setting the initial position and number of stations for the design of a new station network and then use of other methods (see Introduction) to optimize the network;

• (b) use of the proposed methodology as an additional tool/criterion in order to evaluate the position of an existing network which has already been designed by other methods;

• (c) use of the proposed methodology as a basic tool to investigate the changes on grid attributes and cluster formation due to changes in the position and number of stations. This case is very important since it can optimize not only the representativeness of the station network but also the quality of the grids, which is strongly dependent on the data of existing stations.

The application of the proposed statistical techniques captured the spatiotemporal variation of precipitation in Nigeria based on gridded precipitation data. Cross correlation revealed an intense change in monthly latitude, longitude, and elevation dependencies of precipitation. These changes were further investigated using PCA, which provided an excellent visualization of the northward movement effects of IDT on rainfall patterns, identifying also the period of its occurrence (June–September with a peak in August). The cluster analysis segmented the country of Nigeria into territories based on different precipitation patterns, which were further used to evaluate the representativeness of the existing network of rain-gauges to describe the spatiotemporal variation of precipitation. The cluster analysis revealed 47 clusters of which 22 do not have a rain-gauge, while eight clusters have more than one rain-gauge. Thus, more rain-gauges and a better distribution are required to describe the spatiotemporal variability of P in Nigeria. The proposed approach for the evaluation of rain-gauge networks (or meteorological stations in general) is mainly proposed as a complementary methodology and not as a complete method to substitute for other techniques since it is based on a coarse temporal resolution (mean monthly step), which is a limiting factor for capturing climatic peculiarities that may appear at smaller time steps. The combination of the specific methods can be implemented for all climatic variables (e.g., temperature, solar radiation, wind speed, etc.) and can support not only the evaluation of meteorological station networks but also other applications related to decision support systems for various management initiatives which are dependent on climate (e.g., water management in agriculture, flood risk assessment, etc.).