The present study spanned 64 municipalities corresponding to 9 Amazonian states from Brazil. The geographical zone covered approximately from 5.03 ºN to 16.45 ºS and 42.83 to 70.93 ºW. Rainfall data refer to historical monthly mean and annual records collected from 218 rain gauges (Figure 1).

The records ranged from the last 16 years for Rondonopolis (Mato Grosso state) to the last 106 years for Manaus (Amazonas state). The last date for data collection was December 2015. Three timescales were considered: individual months with minimum (September) and maximum (March) rainfall, January-June (rainy season), July-December (dry season) and annual rainfall. Note that January-June and July-December overlap winter and summer periods. However, our main goal was to establish a 6-month timescale separating wet and dry seasons.

We performed an analysis on the potential relationship between the deforestation index and the annual rainfall including all the Amazon states. Most investigations try to search for relationships between deforested areas (km²) and rainfall. From our point of view that approach could be fictitious in this case as different states have different spatial extensions. We used instead the percentage of forest loss which relates forest loss (ha) to the total forest area (ha) for a particular site or state. The forest loss percentages are available at http://www.mongabay.com for 2001-2012 period (Butler, 2014). With that proposal, we averaged the available historical annual rainfall record for each studied Amazonian state and correlated them with the corresponding forest loss index.

Geostatistics incorporates the spatial coordinates of observations $\left(x,y\right)$ in data processing. It is fundamental for the modeling of spatial patterns, prediction at unsampled points and estimation of the uncertainty inherent to those predictions (Goovaerts, 1998). The semivariance function, $\gamma \left(h\right)$ , is one of the key components of geostatistics as it summarizes the spatial variation of the studied variable.

Where $\gamma \left(h\right)$ is the value of the experimental semivariance at each separation distance (lag), h, N(h) is the number of data pairs separated by the lag $h$ , $z\left(x_i\right)$ and $z\left(x_i+h\right)$ are the rainfall values at location $x_i$ and $x_i+h$ , respectively. A plot of ((h) versus $h$ defines the experimental semivariogram. Different models are available for fitting the experimental semivariogram (see for instance McBratney & Webster, 1986). We explored three transitive models (spherical, exponential and Gaussian models).

Equation (2) represents the spherical isotropic model. In this case C₀ is the nugget variance or nugget effect ( $C_0\ge 0$ ), $C$ is the structural variance ( $C\ge C_0$ ), $C_0+C$ is the upper limit of the semivariogram (sill) and A₀ is the distance at which semivariogram reaches a constant semivariance value (correlation range parameter). The Spatial Dependence Degree Index ( $SDDI$ ) can be computed as the ratio between the structural variance, $C$ , which characterizes the variance accounted for by the spatial dependence and the sill, $C_0+C$ . For example, $SDDI<0.25$ suggests weak spatial dependence, $0.25<SDDI<0.75$ corresponds to moderate spatial dependence while $SDDI>0.75$ indicates strong spatial dependence. This is a modification to the ratio nugget/sill variance proposed by Cambardella et. al. (1994). The range is the limit of spatial dependence. The spherical model characterizes a moving average of a randomized process (Kuzyakova et. al. 2001).

Equation (3) is the isotropic exponential model. It differs from the spherical model in the rate at which the sill is reached. It represents autoregressive processes of first order (e.g. Markov or Poisson processes) as the autocorrelation function decays exponentially (Wang et. al. 2010).

Equation (4) corresponds to the Gaussian or hyperbolic isotropic model. It is almost similar to the exponential model but the departure from the nugget variance is smoother.

Block kriging is a robust interpolation method for estimating values of the variable (e.g. rainfall values in this case) at ungauged zones (Burguess & Webster, 1980). The block kriging estimator ( $Z_S$ ) of the average rainfall on a given zone ( $S$ ) is built as a linear combination of the available rainfall data:

Where $Z_S^{*}$ is the estimated kriged value of $Z$ in the area $S$ and ${}_i$ are weighing factors such that:

In this case $n$ is the number of rainfall values in the vicinity of $Z_S^{*}$ . The block kriging estimation is unbiased under the condition:

That is, the mathematical expectation of $Z_S^{*}$ is the same as that of the real $Z_S$ value.

The accuracy of rainfall maps was assessed through a cross-validation process (Davis, 1987). The agreement between estimated and measured values was quantified using mean absolute error ( $MAE$ hereafter) (Voltz & Webster, 1990) and goodness-of-prediction statistics ( $G$ statistics hereafter) (Agterberg, 1984). The MAE statistics is in fact a residual sum:

Where $\widehat{z}\left(x_i\right)$ is the predicted rainfall value at location $i$ .

The $G$ statistics is a sort of balance between kriging and sample mean as potential rainfall predictors:

Where $z\left(mean\right)$ is the mean value of all the observations.

Equation (9) defines three possibilities:

It is useful to recall that due to earth curvature the use of geographical coordinates for large scale spatial analysis is inappropriate. Thus, in order to perform the spatial analysis, geographical coordinates were transformed into rectangular $\left(x,y\right)$ coordinates.

Classical statistics (e.g. first and second order statistics and test of normality) and linear regression analysis were conducted using the Statistica^TM Software Package (StatSoft. Inc., 2011). All the geostatistical analyses were performed using the GS+ Geostatistical Software Package (Gamma Design Software, 2001). In particular, we set the minimum lag class distance to 247.25 km. Block kriging estimation was carried out using $4 x 4$ local grids and 16 neighbors within a radius equal to the range of the semivariogram. We selected a total of six map contour levels.