I'm in a discussion with a colleague who believes that kriging should not be used to interpolate a surface from polygon based data (i.e. with kriging applied to centroids), because this effectively translates the spatial uncertainty inherent in any polygon centroid into data error in the kriged surface. I'm sure he's right in principal but cannot help thinking that the uncertainty-error question should not be an absolute argument against the method give that uncertainty is (a) unavoidable and (b) a function of zone size and the degrees to which the data has already been aggregated. I'd be interested to hear views from others, particularly anyone with good experienced of kriging interpolation.
A key question is what "polygon based data" means? Does it mean the spatial average value, is it a count or a classification.
In the case of a spatial average value it is already well known in the geostatistical literature that the kriging equations are easily modified to allow for non-point support data. Next the question if whether you are trying to interpolate the data to point support or to a larger support size. In either case you need to be able to compute the average value of the variogram (an integral) for pairs of points (one in the polygon and one in the larger support size. In practice this is very difficult to program so one would use numerical integration but even that is difficult unless the polygon is actually a rectangle, i.e. a collection of congruent rectangles. In case the data are counts in polygons see a paper by Carol Gotway, J. American Statistical Assn.
If you treat the data values as being associated with the centroid this will introduce unquantifiable errors and will corrupt the kriging variances.
This is somewhat analogous to the application of kriging in mining where the data values are averages over cylindrical cores (typically 12 inches in diameter and 6-8 long, so-called "diamond drill cores"). The objective is to estimate average values for "blocks", typically 40 ft high, 100 ft by 100 ft.
When you partition a region into polygons, the variability "within" the region is split between the total variabilities within the respective polygons plus the variabilities of the average values of the polygons within the region. Similar to what happens in ANOVA, in geostatistics this is known as "Krige's Theorem"