Harmonization of design-based mapping for spatial populations

The mapping of a survey variable throughout a continuum or for finite populations of units is usually performed from a model-dependent perspective. Nevertheless, when a sample of locations/units is selected by a probabilistic sampling scheme, the complex task of modelling can be avoided by using the inverse distance weighting interpolator and deriving the properties of maps in a design-based perspective. Conditions ensuring consistency of maps can be derived mainly based on some obvious assumptions about the pattern of the survey variable throughout the study region as well from the feature of the sampling scheme adopted to select locations/units. Nevertheless, in a design-based setting the totals of the survey variable for a set of domains partitioning the study region are commonly estimated by traditional estimators such as the Horvitz–Thompson estimator in the case of finite populations or the Monte-Carlo estimator in the case of continuous populations or by related estimators exploiting the information of auxiliary variables. That necessarily gives rise to different total estimates with respect to those achieved from the resulting maps as the sum of the interpolated values within domains. To obtain non-discrepant results, a harmonization of maps is here suggested, in such a way that the resulting totals arising from maps coincide with those achieved by traditional estimation. The capacity of the harmonization procedure to maintain consistency is argued theoretically and checked by a simulation study performed on some real populations.


Introduction
Spatial phenomena frequently need detailed information especially in natural resources management. Therefore, mapping becomes crucial for the understanding of spatial patterns of an interest attribute. The population to be surveyed can be a continuum, a finite collection of areas portioning a study region, or a finite collection of points spread throughout a study region. In most cases, the available resources render impractical the complete census of these populations. Therefore, the survey variable is recorded only for a subset of locations/areas/points and an estimation criterion is adopted to estimate the values of the survey variable within non recorded locations/areas/points and obtain wall-to-wall maps of the interest variable throughout the whole population.
Usually, mapping is approached in a model-based framework: locations/area/points where the variable is recorded are considered as fixed, while values are assumed to be outcomes of a superpopulation probability model (e.g., Cressie 1993). As an alternative approach, Fattorini et al. (2018aFattorini et al. ( , b, 2019 propose mapping continuous populations, finite populations of areas and finite populations of points in a design-based framework. In this scenario, values are viewed as fixed constants and the probability distribution of any sample statistic is determined from the uncertainty entailed by the probabilistic sampling scheme adopted for selecting locations/areas/points. In their works, the authors highlighted as any design-based mapping is challenging. Indeed, when estimating the values of the survey variable for a single location/area/point, either it has been sampled and therefore there is no need for estimation, or it is unsampled and therefore we do not have any information for performing estimation. Consequently, even in a design-based framework, the use of an assisting model, based on some auxiliary information, seems to be the only way to fill the lack of sample information. As a solution, Fattorini et al. (2018aFattorini et al. ( , b, 2019) considered as assisting model the well-known Tobler's first law of geography, i.e., locations/areas/points close in space tend to have more similar values than those that are far apart (Tobler, 1970). Accordingly, the authors proposed the adoption of the so-called inverse distance weighting interpolator (IDW) (e.g., Henley 1981) and determined the conditions ensuring IDW design-based unbiasedness and consistency.
Nevertheless, it should be acknowledged that frequently, as it occurs for example in forest inventories, estimation of totals and averages of interest attributes for a set of domains partitioning the study region is traditionally performed from a design-based perspective adopting the Horvitz- Thompson (HT) estimator in the case of finite populations or the Monte-Carlo estimator in the case of continuous populations or some modifications able to exploit auxiliary information (e.g., Särndal et al. 1992, chapter 6).
However, total estimates for these domains can be achieved also as totals or integrals of the IDW interpolations within. Obviously, the two resulting sets of estimates will differ. Therefore, the aim of this work is to harmonize IDW maps by rescaling the interpolated values in such a way that totals or integrals of single estimates within domains will match the traditional design-based estimates, avoiding unsuitable discrepancies in the final results, that may be perplexing especially in a report phase. It is worth noting that the harmonization procedure is not introduced to improve the design-based performance of the IDW interpolation, but just to avoid discrepancies of the total estimates achieved from maps with those achieved from traditional methods. Therefore, the capacity of the harmonization procedure to maintain the design-based consistency ensured by the original IDW interpolation has been the main concern. Accordingly, the main target of this paper is to theoretically prove the harmonization consistency and to check it empirically by simulations.
Because averages are achieved as totals divided by the size of the study region in continuous populations or by totals divided by the population size in finite populations of areas or points, and because these quantities must necessarily be known to perform mapping, harmonization with respect to totals or averages are equivalent from a statistical point of view. Moreover, regarding the types of populations considered for harmonization, it should be noticed that for point populations over large study regions (e.g., plants, shrubs, or trees), the list and locations of population units are not available. As consequence, maps of these populations are precluded. The sole cases in which mapping is possible occur for forest stands located on surfaces of limited size (few hectares) in which 3P sampling can be performed. Indeed, in these cases it becomes possible to visit (and then to locate and list) all the trees by a team of experts and to give predictions of the survey variable for each tree, that is a necessary step to perform 3P sampling (Gregoire and Valentine 2008). However, in most of these cases, the principal aim is the estimation of the total timber volume, while mapping is less urgent: this is due to the small size of the stand and also because the predictions from a crew of experts are likely to well depict the distribution of the timber volumes throughout the stand. For these reasons, harmonization in finite populations of points is not considered in the paper.
The paper is organized as it follows: in Sect. 2 some preliminary results on design-based mapping of continuous populations and finite populations of areas are given. Then, the procedures for harmonizing IDW maps are described in Sect. 3 and 4 where the design-based consistency of the harmonized maps is theoretically proven, and a pseudopopulation bootstrap estimator of map precision is proposed. In Sect. 5 a simulation study is performed on a set of real populations to empirically check the capacity of the harmonization procedure in maintaining the consistency of the resulting maps. Concluding remarks are reported in Sect. 6.

Preliminary results on design-based IDW interpolation
Consider a study region A, that is supposed to be a connected and compact set of R 2 , and let f be a bounded measurable function related to the values of a survey variable Y and defined on a subset B & A. Moreover, let k Á k be a norm in R 2 and / : ½0; 1Þ ! R þ be a nonincreasing continuous distance function on ð0; 1Þ, with / 0 ð Þ ¼ 0 and A widely applied class of distance functions satisfying (1) is the class of negative powers of order a, given by (e.g., Gong et al. 2014;Noori et al. 2014;Bȃrbulescu et al. 2021). The choice of (2) is particularly appealing owing to its simplicity and because of the straightforward interpretation of the a parameter. Indeed, as showed in the next subsections, the IDW interpolator is a convex combination of the sample values, with weights proportional to the values of the distance function /. Therefore, negative powers of distances obviously give less weight to sample values further away from the location where interpolation is performed with a that plays the role of the smoothing parameter, i.e., for smaller values of a, the interpolated map becomes smoother, till, in the limiting case of a ¼ 0, the map is constant as the interpolated values are all equal to the sample mean. On the other hand, when a becomes larger and larger, the estimated map becomes rougher and rougher till, for a approaching infinity, the IDW interpolator reduces to the well-known nearest neighbour (NN) interpolator (Fattorini et al. 2021).
Depending on the features of the spatial populations, the IDW mapping of the interest attribute is performed in the following two settings.

Continuous populations
B coincides with A and f ðpÞ is the value or the density of the survey variable Y at p 2 B. Therefore, mapping necessitates the knowledge of f ðpÞ for each p 2 B. To this purpose, let P 1 ; . . .; P n be n random variables with values in B that represent the n locations selected from B by means of a probabilistic sampling scheme. Then, in accordance with Fattorini et al. (2018a), if the existence of the continuous probability density function of ðP 1 ; . . .; P n Þ is assumed, for a fixed a, the IDW interpolator of f ðpÞ is almost certainly equal to with weights that under the distance function (2) turn out to be w i;a p ð Þ ¼ kP i À pk Àa P n h¼1 kP h À pk Àa ; i ¼ 1; . . .; n: For a ! 1 the interpolator (3) reduces to the NN interpolator, that is almost certainly equal to where P NNðpÞ ¼ argmin i¼1;...;n kP i À pk. In this case the resulting map is a surface piecewise constant on the Voronoi cells around the sampled locations (Fattorini et al. 2021).

Finite populations of areas
A is partitioned into a finite population U of N spatial units a 1 ; . . .; a N of extents k 1 ; . . .; k N . In this case B is the set of centroids b 1 ; . . .; b N of the areas and y j is the amount of the survey variable Y within a j . Therefore, mapping necessitates the knowledge of y j for each j 2 U. However, since the area extents are usually known, knowledge of the y j s is equivalent to the knowledge of densities f j ¼ f b j À Á ¼ y j =k j for each j 2 U. Accordingly, if S denotes the set of labels identifying the selected areas, then in accordance with Fattorini et al. (2018b), for a fixed a, the IDW interpolator of f j is given by where Z j is the random variable equal to 1 if j 2 S and equal to 0 otherwise, with weights that under the distance function (2) turn out to be where is the set of centroids in the sample that are nearest to b j . Indeed, contrary to the case of continuous populations, in the case of finite populations of areas the nearest neighbours of an area may be more than 1 as, for example, in the case of populations of regular polygons such as pixels (Fattorini et al. 2021).
Once the b f j;a s are achieved, the interpolation of the y j s is simply given by b y j;a ¼ k j b f j;a for j 2 U. However, the interpolation of densities is more suitable for working in the asymptotic scenario considered by Fattorini et al. (2018bFattorini et al. ( , 2021 in which the k j s decrease and then the y j s approach zero.

Consistency conditions
The asymptotic properties of the IDW interpolators (3) and (5) are derived in Fattorini et al. (2018a, b), respectively, while those of their NN counterparts (4) and (6), achieved for a ! 1, are derived in Fattorini et al. (2021). Without going into the technical details provided in those papers, for both the population settings, the design-based unbiasedness and consistency of IDW and NN interpolators concern: i) a sort of smoothness of the function f onto the study region with jumps that occur in sub-sets of measure 0; ii) when dealing with finite populations of areas, the regularities in the area shapes; iii) the capacity of the sampling design to asymptotically achieve a spatial balance of the selected locations/areas; iv) the use of negative power distance functions with a [ 2.
It is worth noting that the four conditions seem to match most real situations encountered in environmental surveys. Indeed, the smoothness assumption i) is very common and it is at the basis of most interpolation techniques (e.g., Cressie 1993, Sect. 3.1). Moreover, assumption i) is also reasonably valid in many natural scenarios where the density of an attribute changes smoothly throughout space (continuity) and when it changes abruptly, that usually occurs along borders delineating variations in the characteristics of the study region (e.g., forest-meadows). Therefore, borders may be realistically approximated by curves well approaching the theoretical condition of discontinuity over a region of zero measure. Regarding assumption ii) concerning the regularity of the shape of areas, it is ensured by the fact that in most cases, especially in forest and vegetation surveys, areas are regular polygons. Finally, iii) and iv) actually do not constitute assumptions because both sampling schemes and distance functions are chosen by the user. In particular, as emphasized in Sect. 3, the asymptotic spatial balance required by iii) is ensured under the schemes usually applied in environmental surveys.

Data-driven choice of a
Because any value of a [ 2, including a ¼ 1, ensures the design-based consistency of the IDW interpolation, Fattorini et al. (submitted) propose to choose a, that in this framework plays the important role of a smoothing parameter, by means of a data driven procedure. In particular, the authors propose the use of the leave-one-out cross validation (LOOCV) that constitutes an intuitive and widely applied technique in spatial interpolation (e.g., Giraldo et al. 2011;Ignaccolo et al. 2014;Montanari and Cicchitelli 2014).
LOOCV consists in removing one location/area at a time from the sampled ones, interpolating the value or density of the survey variable at the removed location/area using all other locations/areas in the sample and then repeating this process for each location/area in the sample. The interpolated values at each sample location/area are then compared with the actual values minimizing the sum of squared differences.
Accordingly, in the case of continuous populations, a is selected to minimize where b f Ài;a P i ð Þ is the IDW interpolator of f P i ð Þ achieved by means of (3) or (4) from the sample of n À 1 locations obtained deleting the sample location P i .
Similarly, in the case of finite populations of areas, a is selected to minimize where b f Ài;a is the IDW interpolator of f i achieved by means of (5) or (6) from the sample of n À 1 areas obtained deleting the area i 2 S. In the following, the IDW interpolators with a chosen by means of LOOCV are denoted by b f b a ðpÞ for each p 2 B in the case of continuous populations or by b f j;b a for each j 2 U in the case of finite populations of areas and are referred to as data-driven IDW (DD-IDW) interpolators. Because b a is chosen from sample data instead of being fixed in advance, it is a random variable whose presence may, in principle, increase the variability of the DD-IDW interpolators and may preclude their consistency. Fattorini et al. (submitted) broad the consistency results achieved for IDW interpolators proving the consistency of the DD-IDW under the same asymptotic scenarios.

Bootstrap estimation of precision
Regarding the estimation of the precision of the DD-IDW interpolators, Fattorini et al. (submitted) propose the use of a pseudo-population bootstrap (see e.g., Quatemberg 2015) based on adopting the DD-IDW map achieved from the sample as the pseudo-population from which M bootstrap samples are re-sampled using the same scheme adopted to select the original sample and then achieving as many DD-IDW bootstrap maps from these samples.
Accordingly, in the case of continuous populations, let the sample values f ðP 1 Þ; . . .; f ðP n Þ. Then, for each p 2 B, the pseudo-population bootstrap estimator of the root mean squared error (RMSE) of b f b a ðpÞ is given by where P Ã 1;m ; . . .; P Ã n;m are the locations selected in the mth bootstrap resampling by means of the scheme adopted to select the original sample be the DD-IDW map based on the sample densities f j ; j 2 S. Then, for each j 2 U, the pseudo- where S Ã m is the sample of areas selected in the m-th bootstrap resampling by means of the scheme adopted to select the original sample, b f j;b a ; j 2 S Ã m are the densities in these areas derived from the estimated map Fattorini et al. (submitted) prove that for large sample sizes and for a sufficiently large M, the bootstrap estimators (9) and (10) tend to be conservative with the ratio of the expectation of the bootstrap RMSE to the true value that is bounded by ffiffiffiffiffi 10 p . Even if this result may induce to suspect a large overestimation of the RMSEs, such bound should be viewed as a threshold limiting possible overestimation.

Harmonization with the overall population total
In environmental and ecological studies, the total of a survey variable for the whole study region frequently covers a role of great interest. As typical examples, estimation of the total amount of a pollutant in a lake is crucial for achieving information on environmental changes in the surrounding zones (e.g., Greaver et al. 2016), estimation of the total erosion extent in a region is fundamental for the management of cultivations (e.g., Kelley, 1990), estimation of the total carbon storage in a forest is essential for determining sequestration capacity (e.g., FAO 2010, Chapter 2). Once locations or areas are selected by means of a probabilistic sampling scheme, the population total is commonly estimated by means of the HT criterion or by related criteria able to exploit the presence of auxiliary information. These criteria constitute consolidated, widely applied strategies as they were unbiased or approximately unbiased with design-based variances with known analytic expressions or approximations. Moreover, variances can be estimated on the basis of those expressions/approximations (Gregoire and Valentine, 2008).
However, total estimates naturally arise also from the resulting maps as the integral, in the continuous case, or the sum, in the discrete cases, of the interpolated values, thus achieving total estimates that invariably differ from those achieved by traditional, HT-based techniques. To eliminate these unsuitable discrepancies, the matching of the two total estimates can be performed by rescaling the DD-IDW interpolations of the interest attribute for each location/ area. Accordingly, for continuous populations, the rescaled DD-IDW map is given bỹ while, for finite populations of areas, the rescaled DD-IDW map is given bỹ respectively, where, in both cases, b T denotes the total estimate obtained by means of a commonly adopted HTbased technique and b T b a denotes the total estimate obtained from the resulting DD-IDW map.

Harmonization for continuous populations
In the case of continuous populations, the spatial total can be expressed as (see e.g., Stevens 1997; Gregoire and Valentine 2008, Chapter 10), while its design-based estimation b T can be obtained extending the HT estimator to the continuous case where pðP i Þ denotes the strictly positive inclusion density function on B at the sample locations P i ði ¼ 1; . . .; nÞ (see e.g., Cordy, 1993).
Because of the integral representation (13), the estimation of T may be approached as a Monte Carlo integration. Interestingly, the most common Monte Carlo integration methods such as crude Monte Carlo integration, modified Monte Carlo integration and random grid Monte Carlo integration are equivalent to Uniform Random Sampling (URS), Tessellation Stratified Sampling (TSS) and Systematic Grid Sampling (SGS), respectively, that constitute the most common sampling schemes adopted in environmental surveys for continuous populations (e.g., Barabesi 2003). In these cases, the estimator (14) reduces to Consistency of (15) under URS, TSS and SGS has been proven by supposing a sequence of designs, each of them selecting an increasing number of locations. In particular, the relative efficiency of TSS with respect to URS has been proven for finite samples, also proving that efficiency approaches infinity as the sample size increases because TSS variance goes to 0 more quickly than under URS (Barabesi and Franceschi 2011;Barabesi et al. 2012Barabesi et al. , 2015. Moreover, consistency of (15) under SGS has been proven by Fattorini et al. (2020).
On the other hand, from the resulting DD-IDW map, the total estimate turns out to be For a fixed and under the same asymptotic scenario and the same schemes (URS, TSS and SGS) that ensure the design-based consistency of (15), Fattorini et al. (2018a) prove the design-based consistency of the IDW estimator b f a ðpÞ that, in turn,-owing to the Lebesgue dominated convergence Theorem -entails the consistency of the Þdp to the true total T. Therefore, stated the consistency of the DD-IDW interpolator b f b a ðpÞ under the same asymptotic scenario, that also entails the consistency of (16).
Joining the two consistency results for the estimators (15) and (16)

Harmonization for finite populations of areas
In the case of finite populations of areas, the population total can be expressed as while its design-based estimation can be obtained by means of the HT estimator where p j denotes the first-order inclusion probability of the area j induced by the sampling scheme adopted to select areas. As the extents of areas partitioning the study region decrease in such a way that their number and sample size increase, Fattorini et al. (2020) derived conditions ensuring design-based consistency of (18). In particular, they proved that consistency conditions are satisfied when Simple Random Sampling Without Replacement (SRSWOR) and One Per Stratum Sampling (OPSS), are adopted. Moreover, they proved consistency also under Systematic Sampling (SYS), that however necessitates further assumptions.
On the other hand, from the resulting DD-IDW map, the total estimate turns out to be For a fixed and the same asymptotic scenario and the same schemes (SRSWOR, OPSS and SYS) that ensure consistency of (18), Fattorini et al. (2018b) prove the design-based consistency of the IDW estimator b f j;a that, in turn-owing to the Result B.1, Appendix B of that paperentails the consistency of b T a ¼ P j2U k j b f j;a to the true total T. Therefore, stated the consistency of the DD-IDW interpolator b f j;b a under the same asymptotic scenario, that also entails the consistency of (19). Joining the two consistency results for the estimators (18) and (19), the rescaling constant b T = b T b a in Eq. (12) obviously converges to 1, in such a way that the harmonized interpolatorf j;b a converges to the DD-IDW interpolator b f j;b a . That proves consistency of the harmonized interpolator (12) for finite populations of areas under SRSWOR, OPSS and SYS.

Bootstrap estimation of precision
Regarding the estimation of the precision of the harmonized interpolators, because harmonization is performed from sample data, rescaling the DD-IDW interpolations by means of the ratio b T = b T b a involves uncertainty that must be accounted in the bootstrap procedure. Therefore, in the case of continuous population, the bootstrap RMSE estimator (9) is changed intõ are the total estimates (18) and (19) achieved from the bootstrap sample values b f j;b a ; j 2 S Ã m .

Harmonization by domains
It frequently occurs that total estimates of an interest attribute are also required for D subpopulations partitioning the entire population, usually referred to as domains (e.g., Särndal et al. 1992 Sect. 10.3). For example, a survey region may be divided into D domains by administrative bounds (e.g., regions, counties, municipalities) and we could be interested, besides to the mapping and to the overall total of an attribute, to its sub-totals within the domains.
In the case of continuous populations, denote by B 1 ; . . .; B D the D subsets partitioning B, whose totals are of interest. In this case, estimation of totals within domains can be performed simply defining, for each domain d ¼ 1; . . .; D, the function.
in such a way that the total for the domain d, say T ðdÞ , is obtained as in (13) by substituting f p ð Þ with f ðdÞ p ð Þ. Similarly, the Monte Carlo estimator for the d-th domain, say b T ðdÞ , can be performed as in (15), once again substituting f p ð Þ with f ðdÞ p ð Þ. On the other hand, from the resulting DD-IDW map, the total estimate for each domain d ¼ 1; . . .; D, say b T b aðdÞ , turns out to be as in (16) with the integral extended to B d , instead of the whole B. Therefore, for continuous populations, the rescaled map that ensures the matching of total estimates for each domain as well as the matching of the overall total estimates is given bỹ Consistency of (22) arises, mutatis mutandis, from the same consideration performed in Sect. 3.1.
In the case of finite populations of areas, denote by U 1 ; . . .; U D the D subsets of areas partitioning U, whose totals are of interest. In this case, estimation of totals within domains can be performed simply defining, for each domain d ¼ 1; . . .; D, the values.
in such a way that the total for the domaind, sayT ðdÞ , is obtained as in (17) (18), once again substituting f j with f jðdÞ . On the other hand, from the resulting DD-IDW map, the total estimate for each domaind ¼ 1; . . .; D, say b T b aðdÞ , turns out to be as in (19) with the summand extended toU d , instead of the wholeU. Therefore, for finite populations of areas, the rescaled map ensuring the matching of total estimates for each domain as well as the matching of the overall total estimates is given bỹ Consistency of (23)

Simulation studies
For each population setting, a simulation study is performed to empirically check if and how much the harmonization procedure deteriorates the performance of DD-IDW maps, as well as to check the rate of convergence of the harmonized maps to the original ones.

Populations and sampling
As to continuous populations, the survey region considered for the simulation was a quadrat region of 90,000 ha located in North-Western Tuscany (Central Italy). The population values to be mapped were the precipitations (mm) occurred between 3rd January 2021 and 3rd February 2021, that were artificially achieved by means of an ordinary kriging prediction performed from the values recorded on 32 rain gauge stations that were present in the survey region. The population average was Y ¼ 299 mm for the whole region (see Fig. 1). Moreover, the region was partitioned into D ¼ 2; 4; 8 domains of equal sizes, as depicted in Fig. 2. Precipitation averages (in mm) within the two domains were Y ð1Þ ¼ 306 and Y ð2Þ ¼ 293, precipitation averages within the four domains were Y ð1Þ ¼ 255, Y ð2Þ ¼ 361, Y ð3Þ ¼ 323, and Y ð4Þ ¼ 257, precipitation averages within the eight domains were Y ð1Þ ¼ 219, Sampling was performed by selecting n ¼ 16; 36; 64; 100 locations on the quadrat region by means of URS, TSS and SGS, i.e., selecting n locations completely at random (URS), partitioning the quadrat into n subquadrats of equal size and selecting a location at random within each of them (TSS), or selecting one location at random in one of them and then repeating the location in the others (SGS).
Sampling was performed by selecting samples of n ¼ N=10 areas by means of SRSWOR, OPSS and SYS, i.e., selecting n areas at random without replacement (SRSWOR), partitioning the populations into n blocks of 2 Â 5 contiguous areas and selecting an area at random within each of them (OPSS), or selecting one area at random in one block and then repeating it in the others (SYS).

Simulation
For each combination of population, sampling scheme and sample size, sampling was replicated R ¼ 10; 000 times. At each simulation run, the DD-IDW map was obtained by  Fig. 3 Maps of the growing stock volumes in a rectangular region located in Calabria (Southern Italy) partitioned into three populations of 250, 1,000 and 4,000 areas means of the LOOCV selection of a, that was performed by minimizing (7) or (8), for the continuous population or finite populations of areas, respectively, starting with a ¼ 3 and increasing a by one. Since for quite large values of a the IDW interpolator is practically indistinguishable from the NN interpolator, if the minimum of (7) or (8) was reached for a ¼ 21, the NN interpolator was used. Moreover, the bootstrap RMSE estimators (9) or (10) were adopted in the case of the continuous population and finite populations of areas, respectively, computed by means of M ¼ 1000 bootstrap samples. Then, at each simulation run, harmonization was performed by estimating the population totals by means of the traditional estimators (15) or (18) in the case of the continuous populations or finite populations of areas, respectively, by estimating the totals from the resulting maps by means of estimators (16) or (19) in accordance with the two types of populations, and then harmonizing the DD-IDW by rescaling them by means of Eq. (11) for the continuous population and by means of Eq. (12) for finite populations of areas. Finally, domain partitions were considered for each population as depicted in Figs. 2 and 4 and harmonization was performed for each of the D ¼ 2; 4; 8 domains as described in Sect. 4. Regarding the estimation of precision for the harmonized maps, the bootstrap RMSE estimators (20) or (21) were adopted for the continuous population or the finite populations of areas, respectively, by means of M ¼ 1000 bootstrap samples.

Performance indicators
In the case of the continuous population, interpolation was performed on a regular grid of 100 Â 100 within the quadrat region of Fig. 1. Let b f r p k;l À Á andf r p k;l À Á be the DD-IDW and the harmonized interpolations of the continuous population computed at the node p k;l of the grid (k; l ¼ 1; . . .; 100) at the r simulation run, and let d rmse Ã r p k;l À Á andrmse Ã r p k;l À Á be the bootstrap RMSE estimates from Eqs. (9) and (20), respectively. Based on the R ¼ 10; 000 runs, the absolute bias (AB) the RMSE and the bootstrap ratio (BORAT) were computed for the DD-IDW interpolation at each node p k;l for k; l ¼ 1; . . .; 100. The same indicators were computed for the harmonized interpolation replacing b f r p k;l À Á withf r p k;l À Á in Eqs. (24) and (25) and replacing d rmse (26). Similarly, in the case of finite populations of areas, let b f j;r andf j;r be the DD-IDW and the harmonized interpolations for the area j in the finite populations of areas at the r simulation run, and let d rmse Ã j;r andrmse Ã j;r be the bootstrap RMSE estimates from Eqs. (10) and (21), respectively. Based on the R ¼ 10; 000 runs, the AB the RMSE and the BORAT were computed for each area j. The same indicators were computed for the harmonized interpolation replacing b f j;r withf j;r in Eqs. (27) and (28)   Finally, let b T r and b T b a;r be the total estimates (overall or within a domain) achieved at the r simulation run by means of the HT or Monte Carlo estimators and from the DD-IDW maps, respectively. The Monte Carlo distribution of the correction factors was considered to evidence the level of matching between the DD-IDW and the harmonized maps.

Results
Tables SI1 and SI6 in the Online Resource file contain the minima, the averages, and the maxima of ABs, RMSEs, and BORATs arising from the DD-IDW interpolation for each population, sampling scheme and sample size.  Tables SI2 and SI7 for the harmonized interpolation with respect to the overall total. Figures showing the spatial pattern of these indicators are not reported because they resulted quite similar to those achieved using the DD-IDW interpolation. Moreover, for each population, sampling scheme and sample size, Tables SI3 and SI8 report minima, averages, and maxima of the Monte Carlo distributions of the correction factor (30). Finally, Tables SI4 and SI9 contain the minima, the averages, and the maxima of ABs, RMSEs, and BORATs arising when harmonization is performed with respect to the total estimates within D ¼ 2; 4; 8 domains, while Tables SI5 and SI10 report minima, averages, and maxima of the pooled Monte Carlo distributions of the RxD correction factors (30) in presence of D ¼ 2; 4; 8 domains.
Simulation results show that harmonized and DD-IDW maps are comparable in terms of ABs, RMSEs and BOR-ATs for all the populations and sampling schemes, suggesting that harmonization can be performed without relevant loss in accuracy, precision, and bootstrap performance. Furthermore, in all the situations, the correction factor (30) quickly approaches to 1 as the sample sizes increase. However, the RMSEs of the harmonized maps tend to increase as the number of domains increases. This is an expected result due to the fact that, when the number of domains is large and their sizes decrease, the number of selected units within domains becomes smaller and smaller, thus invariably reducing the precision of the total estimates from which harmonization is performed.

Final remarks
The recent papers by Fattorini (2018aFattorini ( , 2018bFattorini ( , 2019 propose a novel approach for mapping continuous populations or finite populations of areas and points in a design-based framework, avoiding the complex task of modelling the surfaces or the finite populations to be mapped. However, in the design-based approach, an unresolved problem takes rise. Because totals and averages throughout the whole survey region or within domains are traditionally achieved by the HT or related criteria, these estimates invariably differ from the estimates achieved by the resulting maps. For avoiding these unsuitable discrepancies, that would appear awkward especially in a reporting phase, we here propose to rescale the resulting maps in such a way that the total or average estimates arising from maps will match the traditional estimates. We also prove the asymptotic convergence of the two maps, i.e., the convergence to one of the rescaling factors. That is unambiguously confirmed by the simulation study, as the Monte Carlo distributions of the rescaling factors quickly approach one as sizes increase, even for moderately small domains. Therefore, if we simply look at those results, we may argue that harmonization is a useless procedure because DD-IDW and harmonized maps are very similar. However, if we look at the minima and maxima of the rescaling factor, reported in Online Resource file, it is apparent that in some situations, especially when the estimation of domain totals is involved, the rescaling constant may vary from about 0.1 to 3.4. On the other hand, the deterioration of precision of harmonized maps when domain sizes decrease ought to warn against an acritical use of harmonization. Obviously, as the number of domains increases and their extents become small, the number of sample units available to perform estimation within domains becomes small and, in some cases, may be 0. This fact will necessarily introduce an increase in the variance in the estimator based on the HT or related criteria. As we cannot trust in these estimates, harmonization with respect to small domains should be avoided. These issues are well known in literature as small area estimation problems, for which an extensive body of knowledge has emerged recently (see e.g., Rao and Molina, 2015 and references therein).
At the end of this paper it should be pointed out that the harmonization problem could be completely bypassed if a model-based mapping was adopted. Indeed, in presence of spatial autocorrelation and second-order stationarity of the spatial process that is supposed to generate the population under study, the kriging methods not only provide the best linear unbiased mapping but directly provides the best linear predictor of totals within subareas with no necessity of harmonization (e.g. Thompson, 2012, Sect. 20.4). In addition, a long sequence of alternative model-based mapping procedures is available, such as sandwich mapping (Wang et al. 2013) that is suitable in presence of weak spatial autocorrelation and spatial heterogeneity. Indeed, when dealing with some populations, past experience may have established convincingly that certain types of patterns are typical for the survey variable. In these cases, as pointed out by Wang et al. (2012) the patterns can be used to identify a family of stochastic models (superpopulations) that are supposed to generate the population under study for obtaining the most precise possible mapping for a given amount of sampling effort. Obviously, also in design-based approach, in which the minimal sufficient statistics is the unordered set of distinct labelled observations (e.g. Basu 1969), one would like to be able to say what estimator produces the best mapping. However, because in this case the minimal sufficient statistic is not complete (Thompson, 2012, Sect. 9.4), one cannot make statements about one mapping method being best. The lack of completeness of the minimal sufficient statistic in design-based inference has been the main reason for lack of optimality results in this approach. Therefore, one may wonder why to adopt a design-based mapping with the subsequent necessity of performing harmonization if that problem disappears in model-based approaches with also the possibility of achieving the best mapping. The contraposition between model-based and designbased approaches is a deeply debated issue. Drawbacks and merits of the two approaches are well delineated in both statistical literature (e.g., de Gruijter and ter Braak, 1990;Smith 1994Smith , 2001Little 2004;Thompson, 2012) as well as in environmental applications (e.g., Schreuder et al. 1993;Gregoire, 1998;Gregoire and Valentine, 2008;Wang et al. 2013). However, besides the fact that, as emphasized by Särndal et al. (1992) ''Design-based inference is objective, nobody can challenge that the sample was really selected according to the given sampling design. The probability distribution associated with the design is real, not modelled or assumed'', a further advantage of a design-based approach to mapping includes obtaining consistency of maps just on the basis of the design adopted to select fairly representative or balanced samples of locations, requiring only-as a sort of nonparametric approach-mild assumptions about the population under study. In our case, IDW interpolation only requires a sort of smoothness of the surface to be interpolated with discontinuities that occur in sub-sets of measure 0 and, when dealing with finite populations of areas, regularities in the area shapes. Therefore, our IDW mapping is applicable to all the populations sharing these features. Finally, a very practical motivation is relevant in favour of a design-based approach. As stated before, in environmental and forest surveys estimation of totals and averages avoids model-based procedures and is traditionally performed from a design-based perspective exploiting well experimented sampling schemes such as systematic grid sampling and tessellation stratified sampling (see e.g. Tomppo et al. 2010). Therefore, it would been logically inconsistent to adopt a design-based inference to perform the estimation of totals while adopting a model-based inference for mapping. That ultimately motivates the use of design-based mapping with the subsequent, practical necessity of harmonization.
Funding The authors have not disclosed any funding.

Declarations
Conflict of interest All authors declare that they have no conflict of interest.
Consent to participate All authors have consented to participate.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/.