Philippe Cantet
email: philippe.cantet@irstea.fr

Mapping the mean monthly precipitation of a small island using kriging with external drifts

Mapping the mean monthly precipitation of a small island using kriging with external drifts P.

The present study focuses on the performance of spatial techniques to map the annual and monthly mean precipitation in a small island (1128 km 2 ) named Martinique, located in the Lesser Antilles (see Fig. 1(a)). It is well known that, in this region, precipitation is strongly inuenced by orography and the proximity of land [START_REF] Smith | Orographic precipitation in the tropics: The dominica experiment[END_REF]) leading to great spatial variation. Using predictors for mapping precipitation therefore seems to be appropriate.

With a digital elevation model (DEM), 17 covariates which are likely to explain precipitation were built. Most of the above-mentioned methods were tested via a cross-validation procedure on 35 meteorological stations. For the regression methods, predictors were chosen by established techniques whereas a new approach is proposed here to select external drifts in kriging based on a stepwise model selection by the Akaike Information Criterion (AIC). Because of the limited number of stations, local techniques presented in [START_REF] Sun | Mapping the mean annual precipitation of china using local interpolation techniques[END_REF][START_REF] Xie | Performance evaluation of interpolation methods for incorporating rain gauge measurements into nexrad precipitation data: a case study in the upper guadalupe river basin[END_REF] were not tested.

The paper is sectioned as follows. A description of the data used (observation and predictors) is provided in Section 2. This is followed in Sections 3.1 and 3.2 by the interpolation techniques used in the study. The procedures to evaluate the performance of these methods are mentioned in Section 3.3. Then, selection the predictors is described in Section 3.4. Section 4 presents the performance of the dierent spatial methods in which two cross validation procedures were performed. Finally, results are discussed in Section 5, followed by a summary in Section 6.

Data

Generally, a mapping method requires two types of data. The rst one, which is essential, corresponds to observation of the phenomenon we want to map. The second, which may or may not essential depending on the method used, is some variables, such as elevation, that are known over the whole territory and which can explain the phenomenon; these variables are called predictors in the following.

Observation

The precipitation data used in this study were provided by the French meteorological Institute (Météo-France). This dataset contains monthly precipitation amounts from 1991 -2010. The rain gauge network consists of 37 rain gauges that are distributed across the whole island (see Fig 1(b)) for a surface of 1128 km 2 . Note that for January and December, 2 stations recorded no data during several days and were therefore removed from the database of mean monthly precipitation. Only weather stations with complete records were included.

The main statistical values of precipitation are shown in Tab. 1 to summarize the spatial and temporal variability of mean precipitation. The annual cycle comprises two seasons : the dry period between February and April and the wet season between June and November in which monthly precipitations are at least twice as high as in the dry season.

The monthly precipitation can also vary considerably spatially (e.g. ratio > 6 in January). This variability can be explained by the complex topography of the island of Martinique (See 1(c)). It should be mentioned that the locations of the stations are not ideal: they are lacking in the mountainous area and are unevenly distributed on some coasts. 1 Descriptive statistics of the monthly and annual precipitation (mm). The bottom line indicates, for each month, the number of stations for which the monthly precipitation is known. The last two columns give descriptive statistics of elevation for the 37 stations and for the Shuttle Radar Topography Model (SRTM) data (134142 points).

Building candidate predictors : physiographic variables

Building candidate predictors depends on the phenomenon to be mapped and is a fundamental step.

It is well known that the amount of rain that falls in the West Indies varies immensely from beach to mountain top, so the distance to the coast seems to be a good predictor for precipitation. Three distances were taken into account and are illustrated in Figure 2(g): Besides, Martinique, like most of the islands of the Lesser Antilles, is characterized by a complex topography which tends to enhance the rainfall due to synoptic disturbances by orographic eects [START_REF] Smith | Orographic precipitation in the tropics: The dominica experiment[END_REF]. This means that variables representing the topography should be a good way to predict precipitation.

A high resolution Digital Elevation Model (DEM) is needed to capture the complex topography of the island, since with a too coarse DEM (for example a 1 × 1km 2 resolution), mountains named locally pitons are smoothed, and, hence the elevation of stations provided by the DEM and by metadata can dier signicantly. That is why the SRTM data (resolution ∼ 90m) which is available at http://www2.jpl.nasa.gov/srtm/ was used.

In [START_REF] Daly | A knowledge-based approach to the statistical mapping of climate[END_REF], smoothed altitude was also tested in the precipitation model selection process. Three predictors are introduced and illustrated in Figure 2(a):

z90m : elevation at 90m resolution (SRTM data), zsmooth1 : Average elevation in a 1km neighborhood, zsmooth5 : Average elevation in a 5km neighborhood. [START_REF] Agnew | Gis-based construction of baseline climatologies for the mediterranean using terrain variables[END_REF] emphasized that the maximum elevation in a wedge is very useful in areas with dominant wind directions, because the leeward slopes are drier and warmer than windward slopes. In the present study, the North and South directions were not used since the main perturbations occurring in Martinique come from the Atlantic Ocean (East of the island). Six other predictors were added : zxw1, zxw5, zxw10 : Maximum elevation in the West direction in a 1, 5, 10 km neighborhood (see Fig. This set of 15 physiographic variables was generated from the DEM coming from the SRTM data and explains, a priori, spatial variability in the climate data. These variables can be easily generated with a geographical information system. In our case, the R functions focal and terrain of the raster package [START_REF] Hijmans | raster: Geographic data analysis and modeling[END_REF] were used.

Geographical positions (coordinates X and Y in the UTM 20N spatial reference) were also considered as predictors. 3 Interpolation techniques used in this study First, the interpolation techniques used in this study are briey presented and are separated into two groups: techniques for which predictors (see Sec. 2.2) were used or not. In the following, the notations used are (here, the variable is the mean monthly precipitation):

n the number of sites where the variable is known (observed), s i : a site with observation, i = 1..n, z(s i ) : the amount of the variable at site s i , s 0 : the location where we want to predict the amount of the variable, ẑ(s 0 ) : the value predicted by the tested interpolation method at site s 0 , |s i -s j | the Euclidean distance between s i and s j , X j (s) : the amount of the predictor j at location s.

Without predictors

To apply the following methods, only observation and location of the variable are needed.

1. Inverse distance weighting of order d:

ẑ(s 0 ) = n i=1 λ i .z(s i ) with λ i = 1/|s i -s 0 | d n i=1 1/|s i -s 0 | d (1)
In our case, d = 2 and this method is noted IDW 2. Kriging (simple, ordinary and universal):

Let Z * (.) be the linear regression estimator (Z * (s 0 ) = ẑ(s 0 )) and Z(.) a random function where Z(s i ) = z(s i ).

Then, the value, ẑ(s 0 ), can be written :

ẑ(s 0 ) = µ(s 0 ) + n i=1 ω i (s 0 ). z(s i ) -µ(s i ) (2)
where µ(s 0 ) and µ(s i ) are the expected values of Z * (s 0 ) and Z(s i ), respectively. The kriging weights, ω i (.) must be determined to minimize the estimation variance, Var[Z * (s 0 ) -Z(s 0 )], while ensuring the unbiasedness of the estimator, E[Z * (s 0 ) -Z(s 0 )] = 0. Under some assumptions, kriging is regarded as the best linear unbiased estimator. To simplify, kriging can be based on a linear model for which the residuals are assumed to be spatially dependent.

First, the trend, µ(.), of the random function Z(.) must be specied among 3 types :

µ(s) = m a known constant all over the studied area; this is simple kriging denoted SK µ(s) = µ a unknown constant which is considered to uctuate locally maintaining stationarity within the local neighborhood; this is the ordinary kriging denoted OK µ(s) smoothly varies within each local neighborhood and is modelled as a linear combination of functions of the spatial coordinates; this is the universal kriging denoted UK Then, a theoretical variogram is tted on the empirical variogram determined by the sample values to determine the weights, ω i (.). They are generated when the corresponding system of linear equations, depending on the type of kriging, is solved.

These two methods are exact interpolators, that is to say, ẑ(s i ) = z(s i ) ∀i = 1..n.

Using predictors in interpolation

Contrary to the previous methods (IDW, SK, OK, UK), these methods use auxiliary information to explain the variable we want to map. For example, as temperature commonly depends on the elevation, the goal is then to take into account elevation, called a predictor, to estimate the temperature. Ideally, predictors must be known all over the domain. In this way, 4 methods which can use auxiliary information as predictors are proposed :

Let p be the number of predictors selected.

1. Multiple linear regression (MLR):

ẑ(s) = a 0 + p j=1 a j .X j (s) (3)
where the coecients, a j are determined by the sample values and the least square method. Note that this method imposes n > p. It is recommended to have n > p + 40.

Principal component regression (PCR):

The rst step is to perform a principal component analysis (PCA) on the p predictors. The q chosen eigenvectors P C j are then used as predictors in a MLR :

ẑ(s) = b 0 + q j=1 b j .P C j (s) (4)
where the coecients, b j are determined by the sample values and the least square method. This method is often used instead of MLR when n < p.

Partial least square regression (PLSR):

This method bears some relation to principal component regression, since the method consists in maximizing the variance of predictors (X i ) and the correlation between predictors and the variable of interest.

External drift kriging (EDK):

This method is another type of kriging for which the trend µ(.) is a linear combination of functions of predictors or external drifts. Note that UK is a particular case of EDK where the external drifts are only the spatial coordinates.

The three types of regression (MLR, PCR, PLS) are inexact interpolators, since the values ẑ(s i ) predicted by these models dier from observations z(s i ). The errors ẑ(s i ) -z(s i ), called residuals, follow a gaussian distribution centered at 0. A simple kriging (SK), with µ = 0, can be performed on residuals leading to an exact interpolator. In most cases, kriging residuals signicantly improves the results of the method [START_REF] Nalder | Spatial interpolation of climatic normals: test of a new method in the canadian boreal forest[END_REF][START_REF] Bessa | A method for mapping monthly solar irradiation over complex areas of topography: Réunion island's case study[END_REF]) and is usually called regression-kriging (RK). In our case, three types of regression were used and the fact of kriging residuals will be noted MLRK, PCRK and PLSK.

Remark 1: Although normality is not a prerequisite for kriging, it is a desirable property. Kriging will only generate the best absolute estimate if the random function ts a normal distribution. In this way, applying a mathematical function such as the natural logarithm or power function can be performed on data.

Remark 2: In theory, the predictors for MLR or EDK must be independent. Nevertheless, in practice, this assumption is not really veried. This is the case in the present study where an exhaustive list of predictors was preferred.

Remark3: EDK and MLRK seem to be similar but they lead to dierent results [START_REF] Hengl | Comparison of kriging with external drift and regression-kriging[END_REF]. Contrary to EDK, the parameters of a MLR are estimated by the least square method in which spatial dependence is not taken into account. With EDK, equations are solved at once whereas MLRK explicitly separates trend estimation from spatial prediction of residuals. Note that, for MLRK, there is no danger of instability as is the case with the EDK system.

Besides, in theory, regression requires independent residuals but kriging relies on dependent residuals. For this reason, generalized linear models can be an alternative but this method was not used here.

Choosing the spatial prediction method

Several interpolation methods were presented in Sections 3.1 and 3.2 with the pros and the cons related to each option. From the same data, these methods can provide signicantly dierent maps of the variable of interest (See Fig 3). Choosing the best method is therefore a key point in order to propose a map close to reality.

To achieve this, cross-validation is widely used. The goal is to estimate the expected level of t of a model to a data set that is independent of the data that were used to train the model. It consists in splitting the dataset into training and validation data. For each split, the model is tted to the training data, and predictive accuracy is assessed using the validation data. There are dierent types of cross validation. In our case, only two were used :

1. The Leave-One-Out (LOO) procedure involves using only one observation as the validation set and the n -1 remaining observations as the training set. It requires learning and validating only N split = n times. Then the n ẑ(s i ) predicted values are compared to the n ẑ(s i ) observed values. 2. The repeated random sub-sampling validation involves using randomly p > 1 observations among n as the validation and the n -p remaining observations as the training set. In this study, the number of random splits was xed at N split = 100. Then the (N split × p) ẑ(s i ) predicted values are compared to the (N split × p) z(s i ) observed values.

The best model is the one which minimizes the error between predicted and observed values. To determine this, three criteria were used and are dened below:

Let be z = 1 n n i=1 z(s i ), and, ẑ (s i ) = ẑ(s i ) -z, and, z (s i ) = z(s i ) -z
The Root Mean Square Error (RM SE) is expected to be close to 0.

RM SE = 1 n n i=1 (ẑ(s i ) -z(s i )) 2 (5)
The Nash and Sutclie coecient of eciency (EF ≤ 1) is expected to be close to 1 [START_REF] Nash | River ow forecasting through conceptual models part ia discussion of principles[END_REF].

EF = 1 - n i=1 (ẑ(s i ) -z(s i )) 2 n i=1 (z -z(s i )) 2 (6)
The Willmott index (0 ≤ D ≤ 1) is expected to be close to 1 [START_REF] Willmott | A rened index of model performance[END_REF].

D = 1 - n i=1 (ẑ(s i ) -z(s i )) 2 n i=1 ( ẑ (s i ) -|z (s i )|) 2 (
AIC = 2p -2 ln(L) ( 8
)
where p is the number of predictors and L the likelihood of the model. Given a set of candidate models for the data, the preferred model is the one with the minimum AIC value. In the case n/p < 40 (where n denotes the sample size), it is strongly recommended to use the corrected AICc:

AICc = AIC + 2p(p + 1) n -p -1 (9)
1 Here, a model is dened by predictors

The AIC can be linked with the RM SE criterion [START_REF] Aertsen | Performance of modelling techniques for the prediction of forest site index: a case study for pine and cedar in the taurus mountains, turkey[END_REF]):

AIC = n ln(RM SE) + 2p (10)
Instead of testing all possible models, stepwise methods are widely used to select the best model. In our case, forward selection was used. This involves starting with no predictor in the model, testing the addition of each predictor using a chosen model comparison criterion, adding the predictor (if any) that improves the model the most, and repeating this process until none improves the model. The comparison criterion diers from spatial prediction methods :

for MLR (2 p dierent models), the function stepAIC from the R package MASS was used [START_REF] Venables | Modern Applied Statistics with S[END_REF].

for PLS and PCR (p dierent models), the exhaustive method was used : the LOO procedure is applied for each model to estimate the RM SE. Then the model which minimized the AICc was considered as the best model.

for EDK (2 p dierent models), forward selection was applied with the AICc as selection criterion (evaluated by the RM SE) in performing the LOO procedure for each tested model. The initial model is equivalent to OK.

Remark: For the MLR, PLS and PCR methods, procedures to select predictors are widely used in the literature. For the EDK, the proposed procedure is, as far as we know, a novelty.

Implementing the methods

All the methods were implemented in the software R (R Core Team 2012). For IDW and kriging, the package gstat [START_REF] Pebesma | Multivariable geostatistics in s: the gstat package[END_REF]) was used.

The package pls [START_REF] Mevik | pls: Partial Least Squares and Principal Component regression[END_REF] provides functions for regression methods such as PCR and PLS.

After predicting the mean annual and monthly precipitation with these methods, a set of map layers in raster format was generated. These layers were based on grids, where each point represents the center of a 90m side square (the resolution of the DEM, see Section 2.2 above). Finally, the raster package [START_REF] Hijmans | raster: Geographic data analysis and modeling[END_REF]) was used to illustrate predicted maps (rasters) and to build predictors.

Theoretical and empirical models in kriging

The Matern model with the Stein representation [START_REF] Stein | Interpolation of spatial data: some theory for kriging[END_REF]) was imposed for the kriging for which the semivariance is

γ(h) = c 0 + c 1 1 - 1 2 ν-1 Γ (ν) 2ν 1/2 h r ν Kν 2ν 1/2 h r ( 11
)
where h is the separation distance, Kν is the Basset function, Γ is the gamma function, r is the range or distance parameter (r > 0) which measures how quickly the correlations decay with distance, and ν is the smoothness parameter (ν > 0). c 0 corresponds to the nugget variance, and c 0 + c 1 to the sill variance. The Matern model was chosen since it has great exibility for modelling the spatial covariance [START_REF] Minasny | The matérn function as a general model for soil variograms[END_REF]. Note, that the exponential model γ(h) = c 0 + c 1 exp(-h/r) is a particular case of the Matern model (ν = 1/2).

The empirical variogram was calculated by the method of moments :

γ(h) = 1 N (h) N (h) i=1 [z(s i ) -z(s i + h)] 2 (12)
where N (h) is the number of pairs of observations separated by distance h.

Then, the theoretical model is automatically tted to the empirical variogram by the weighted least squares by minimizing the objective function :

O(θ) = d N d d 2 [γ(d) -γ θ (d)] 2 (13)
where N d is the number of pairs for lag d (see Table 3), γ the empirical variogram, and, γ θ the Matern model with parameter vector θ.

Note that, when a non-stationary (i.e. non-constant) mean is used as in UK and EDK, both for simulation and prediction purposes, the variogram model dened is that of the residual process, not that of the raw observations. In this study, no nugget eect and no anisotropy 2 has been considered. For the interpolation neighborhood, global prediction was used, meaning that for each prediction all the data values were used. In fact, kriging was performed with the default options of the gstat R package. Table 3 The number of point pairs and the average distance (in meters) of all point pairs considered to estimate variograms (in kriging).

because of the lack of points

Results

All the methods described in Sections 3.1 and 3.2 were performed. Residual kriging improved the performance of the methods MLR,PCR and PLS (not shown here). Kriging with only z90m as external drift was also performed, noted z90mK. A log-transformation was also applied on observations in order to obtain more Gaussian data. In order to simplify the tables and gures, the results of only seven methods are mentioned further: IDW, UK, MLRK,PCRK, PLSK, z90mK and EDK.

Predicting performance via the LOO procedure

Table 4 illustrates the performance of the dierent predicting methods for the annual mean precipitation. First, the gain of the log-transformation depends on the methods and is not really signicant. The methods with no predictors (IDW and UK) show the worst results and using them on such a territory seems to be inappropriate. Using predictors signicantly improves the prediction. Adding only the altitude in kriging (z90mK) provides better results than UK but does not appear sucient to fully explain the precipitation, as adding other predictors greatly improves the predicting values. For example, MLRK in which the model is composed of seven predictors shows better results than kriging with altitude. For PLSK, four eigenvectors suce to minimize the AICc. The regression kriging methods (MLRK, PCRK, PLSK) have very similar good performances (RMSE∼280 and EF∼0.87). EDK, which uses only 2 predictors (see Table 5), seems to be the best method to predict mean annual precipitation (RMSE=223.3, EF=0.91). For information, Figure 4 illustrates the empirical and the theoretical variograms for The LOO procedure was also carried out for the twelve monthly mean precipitations (with no log-transformation)

and the results are shown in Fig. 5 . As for the annual precipitation, EDK seems to be the best method to predict mean monthly precipitation,since for any month, the RMSE is the lowest with this method, whereas IDW and UK are not appropriate. For other two criteria (EF and D), there are only three months in which MLRK seems to predict monthly precipitation better. Nevertheless, EDK generally provides the best results. Note that, with a log-transformation of data, EDK is much better and becomes the best method for any criterion and month.

Predictors chosen for EDK

Table 5 lists the predictors used in EDK for the dierent months. The selected predictors dier from one month to another. The descriptive variable zxe1 is the one most often used to predict mean precipitation (9 out of 13 for the Matern variogram), followed by dist_coast and zsmooth5. Note that z90m is never used; in addition, the choice of the function to t the variogram is introduced only from Section 4.5 below.

The number of predictors varies between one and three leading to a model in which the predictors are not strongly correlated. Using the AICc criterion limits this number (See Fig. 6) while minimizing the RMSE should involve choosing eight predictors for annual precipitation, only two are necessary with the AICc criterion. Note, that for PCRK and PLSK, the optimal number of eigenvectors ranges from 4-8 depending on the month. 3 The residuals, after taking into account the dierent trends, were used instead of the target variable to estimate variograms The empirical (points) and the theoretical (lines) variograms for UK (olive green and triangles pointing upwards), z90mK (in pink and triangles pointing downwards), and EDK (in dark red and crossed squares) for the mean annual precipitation. The empirical variogram was estimated and tted from the residual process, not from the raw observations. q q q q q q q q q q q q q IDW UK MLRK q q q q q q q q q q q q q IDW UK MLRK q q q q q q q q q q q q q IDW UK MLRK PCRK PLSK z90mK EDK Fig. 5 Skill scores dened in 3.3 of the dierent interpolation methods estimated by the LOO-procedure for the mean monthly precipitation.

Note that the complexity of the model (characterized by 2p in Eq. 10) has a dierent weight for the total AIC depending on the observation range (which heavily impacts the RM SE) of the target variable. Nevertheless, the number of predictors is the highest for February, the month in which the RMSE is the lowest. This comes from using the logarithm of the RM SE which reduces the scale eect.

× × × × × × × zsmooth1 × × zsmooth5 × × × × × zxe1 × × × × × × × × × zxe5 × zxw10 × slope × Y × Table 5
External drifts (predictors) chosen to predict the mean monthly mean precipitation with EDK. The symbol × (resp.

) indicates that the predictor is selected when the Matern (resp. exponential) variogram is imposed. For PCRK (resp. PLSK), the number of eigenvectors was xed at ve (resp. four) (the optimum number according to the LOO procedure). For EDK and MLRK, their models (predictors) were xed at the beginning (dened by the LOO procedure). Only the variograms or coecients were re-estimated for each random split.

Figure 7 illustrates the skill scores of the seven methods according to the size of the training set to predict the mean annual precipitation. For each method, the performance decreases with the size of the training set, as expected.

Nevertheless, some methods seem to be more sensitive to the number of available data. For example, the performance of IDW degrades only slightly, whereas MLRK, from a critical value, becomes the worst method. So applying this method which showed good performances with the LOO procedure seems to be questionable, probably because the number of stations is not much greater than the number of predictors. The behavior of PCRK and PLSK is very similar while z90mK seems to be more robust than EDK. EDK appears to be the best method for a training set greater than 20, then PLSK becomes better.

Regionalization of the mean annual and monthly precipitation

Predictors are known on the whole territory at a 90m resolution. So, EDK can be performed in order to map the mean monthly (or annual) precipitation at the same resolution in Martinique. These maps are illustrated in Fig. 8 (for annual, February and September) and show the highest amounts of precipitation are located in the regions located at the highest altitude. These amounts are questionable because of the gap between maximum DEM and station altitude (see Table 1). Indeed, the predicted precipitations are extrapolated by EDK to regions with a higher altitude than the highest altitude station (510m). A historical station, which was located at an elevation of 1010m (resp. 965m) according to metadata (resp. DEM), documented the amount of annual precipitation during ve complete years (1962)(1963)(1964)(1965)(1966).

The observed mean annual precipitation at this high altitude station was 4642mm versus 5015mm from the predicted map. This gap of 10% is relatively small but it cannot validate the extrapolation of precipitation. It only shows that the predicted values at high altitude seem to be of the same order of magnitude as what can expected in reality.

4 n = 35 and the size of the validation set is p with p=1..20. size of validation set q q q q q q q q q q q q q q q q q q q q q IDW UK MLRK size of validation set q q q q q q q q q q q q q q q q q q q q q IDW UK MLRK size of validation set q q q q q q q q q q q q q q q q q q q q q IDW UK MLRK PCRK PLSK z90mK EDK In order to study the sensitivity of the method to this choice, the same procedure (selecting predictors by the proposed method, estimating skill scores, mapping) was performed by imposing the exponential model. For most months (9 out of 13), the selected predictors are the same for both models (see Tab. 5). According to the LOO-procedure, the performances are similar : RM SE = 219.5 with an exponential model versus 223.3 with a Matern model for the mean annual precipitation (and log-transformation) or EF = 0.88 with exponential model versus 0.89 with the Matern model for the mean precipitation in February (and log-transformation). Even if the selected predictors are not the same between the two models, the performances remain similar (RM SE = 23 versus 26 for September). The impact of the model choice on the predicted maps are illustrated in Fig. 9. When the predictors are identical for both models, the gap between the two models does not exceed 10% (see Fig. 9(a) and 9(b)) whereas this gap rises to 20% when the predictors dier between the two models (see Fig

9(c))

.

Discussion

The choice of interpolation method depends on data type, area of interest and the spatial scale used. Here, we focused on the mean monthly and annual precipitation in a small island with a complex topography and great spatial variation.

Several spatial methods which are widely used to interpolate mean climate data were tested via cross-validation. The principal novelty of the study resides in the predictor selection in EDK. In the literature, the external drifts are often chosen according to a a-priori knowledge. Selecting predictors according to the proposed method improves the predictive performance of EDK : EF = 0.91 versus EF=0.78 if elevation is the only predictor for the mean annual precipitation (estimated by the LOO-procedure). The gain is similar for all criteria, not only RMSE which is the criterion used to select predictors. The proposed procedure achieves to an accurate prediction of mean precipitation with only one, two or three covariates. Besides, the chosen predictors are quite stable for each pluviometric regime. For example, the months during the dry season (February to April) are characterized by the similar predictors : zsmooth5 and zxe1 even if the selection procedure for these months was performed independently, whereas in the wet season (June to October), dist_coast is chosen for each month. This indicates that there is a climatological coherence in the choice of predictors. Lastly, EDK seems to outperform regression methods such as MLRK, PCRK or PLSK. Overall, the results for mean annual precipitation are notably good when compared to results reported by other studies [START_REF] Hong | Spatial interpolation of monthly mean climate data for china[END_REF][START_REF] Sun | Mapping the mean annual precipitation of china using local interpolation techniques[END_REF]. Moreover, the proposed method to select predictors, performances, and the resulting maps seems to be little inuenced by the choice of the function to t the variogram. The methods used in this paper have also been tested on Guadeloupe (a very similar island to Martinique) leading to similar results : EDK was the most accurate method according to the LOO-procedure and the chosen predictors were coherent and not numerous.

Clearly, this method seems to be a great opportunity to improve the accuracy of interpolations in such islands.

Applying this method to a larger territory with a large number of observations may not be straightforward: the computation is more time-consuming than with other methods since a LOO-procedure must be performed for each predictor to select the best one. Nevertheless, it has been used in France to map the parameters of a rainfall generator [START_REF] Carreau | Extreme rainfall frequency analysis in the south of france : Comparisons of three regional methods[END_REF] with about 2000 rain gauge stations. Observations were split into four climatic zones [START_REF] Cantet | Using a rainfall stochastic generator to detect trends in extreme rainfall[END_REF]) in which the method was performed independently with a reasonable CPU-time. For regions with a (relatively) complex topography such as Mediterranean and Highland, EDK clearly provides the most accurate prediction according to a repeated random sub-sampling validation (as in Section 4.3 above) whereas, in atter regions, the gain of EDK is less pronounced.

Conclusion

In this paper, several spatial interpolation techniques were applied to map the annual and monthly mean precipitation in a small complex topography territory. In our case, if physiographical variables are considered in spatial interpolation, the prediction performances are improved. The cross validation shows that methods which incorporate predictors (MLRK, PCRK, PLSK, z90mK and EDK) clearly outperform methods which consider only coordinates (IDW and UK). This result can be explained by the variable of interest which is a long time scale precipitation estimated over several years. Regression methods (MLRK, PCRK ans PLSK) seem to produce more accurate estimates than kriging with only elevation as drift. However, according to the LOO-procedure, if the predictors selected by the proposed procedure are incorporated in kriging, then EDK gives the best prediction. Note, that the choice of the function used to t the variogram has no signicant impact on the selected predictors, performances or resulting maps of EDK. Nevertheless, the fact that the proposed method is inuenced by the size of the observed sample (more than PLSK but less than MLRK) means that the precautionary principle must be applied. Further investigation is needed to completely validate this method. An alternative method might be to reinterpret distance in IDW (which is the method that is the least inuenced the sampling) as in [START_REF] Ahrens | Distance in spatial interpolation of daily rain gauge data[END_REF]).

Fig. 1

 1 Fig. 1 Location of Martinique in the West Indies (a). Location of the 37 stations in Martinique (b). Topography of Martinique provided by SRTM data (c).

  2(e)), zxe1, zxe5, zxe10 : Maximum elevation in the East direction in a 1, 5, 10 km neighborhood (see Fig.2(c)),[START_REF] Daly | A knowledge-based approach to the statistical mapping of climate[END_REF] also suggested using some characteristics of the terrain since the relationship between precipitation and elevation can vary from one slope to another, depending on location and orientation : slope : slope evaluated with 4 neighbors (see Fig. 2(b)), curv : curvature evaluated with 4 neighbors (see Fig. 2(d)), aspect : aspect evaluated with 4 neighbors (see Fig. 2(f )),

  Fig.2Candidate predictors used in this study.

  Fig. 3 Mean annual precipitation distribution according to three dierent interpolation methods.

  Fig.4The empirical (points) and the theoretical (lines) variograms for UK (olive green and triangles pointing upwards), z90mK (in pink and triangles pointing downwards), and EDK (in dark red and crossed squares) for the mean annual precipitation. The empirical variogram was estimated and tted from the residual process, not from the raw observations.

  Fig.6Values of RMSE (in blue, scale on the left axis) and AICc (in red, scale on the right axis) during the forward selection procedure for the mean annual precipitation (in Fig.6(a)) and the mean precipitation in November (in Fig.6(b)). OK (ordinary kriging) indicates no predictors.

Fig. 7

 7 Fig. 7 Skill scores dened in Section 3.3 of the dierent interpolation methods according to dierent sizes of training set for the mean annual precipitation.

Fig. 9

 9 Fig. 9 Sensitivity of the choice of variogram on the mean annual and monthly (February and September) precipitation: Delta (in percentage) between the maps provided by EDK with a Exponential variogram and with the Matern variogram.

Table 2

 2 Pearson correlation coecient between predictors and the mean annual precipitation (estimated with 35 values).The Akaike Information Criterion (noted AIC) can be used to select a model 1 . Unlike the criteria presented in Section 3.3, it handles the trade-o between the goodness of t of the model and the complexity of the model (number of predictors) and is written as follows:

	aspect	curv	dist_atlan	dist_caraib	dist_coast	slope	z90m	zsmoth1
	0.11	0.07	0.21	0.23	0.79	0.29	0.83	0.87
	zsmooth5	zxe10	zxe1	zxe5	zxw10	zxw1	zxw5
	0.84	0.52	0.79	0.67	0.81	0.80	0.73

3.4 Predictor selection

In Section 3.2, spatial prediction methods which use predictors were presented. In our case, 17 predictors were selected to explain precipitation (see Section 2.2). The fact of taking a predictor into account or not has a signicant impact on the prediction. The Pearson correlation coecient is widely used to show how a predictor is linked with the variable of interest (see Tab. 2). However, choosing a model with predictors which have the greatest correlation is not suitable because two predictors can be strongly correlated (e.g. z90m and zsmooth1) leading to a redundancy of information which can penalize the mapping. Adding one predictor or another is a critical step in mapping.

Table 4

 4 UK, z90mK, and EDK 3 . Residuals coming from the trend of EDK seem to be more ttable than those of UK or z90mK. Skill scores dened in Section 3.3 of the dierent interpolation methods estimated by the LOO-procedure for annual mean precipitation. The cells in grey indicate the best predicting method according to the score.

		With log-transformation	No transformation
		RMSE	D	EF	RMSE	D	EF
	IDW	543.7	0.69	0.56	532.0	0.69	0.58
	UK	484.0	0.73	0.65	495.5	0.73	0.63
	MLRK	284.4	0.84	0.88	269.3	0.84	0.89
	PCRK	295.0	0.84	0.87	290.3	0.84	0.87
	PLSK	295.2	0.84	0.87	275.0	0.85 0.88
	z90mK	386.8	0.778	0.782	388.7	0.77	0.77
	EDK	223.3	0.85	0.91	235.0	0.85	0.91