CSIS Discussion Paper No. 9
Accuracy of Areal Interpolation:
A Comparison of Alternative Methods
Yukio Sadahiro
FEBRUARY, 1999

Center for Spatial Information Science and Department of Urban Engineering
University of Tokyo
7-3-1, Bunkyo-ku, Tokyo 113-8656, Japan

February 12, 1999
Accuracy of Areal Interpolation: A Comparison of Alternative Methods

Yukio Sadahiro
Center for Spatial Information Science and Department of Urban Engineering
University of Tokyo
7-3-1, Hongo, Bunkyo-ku, Tokyo 113-8656, Japan
Phone: +81-3-3812-2111 (ext. 6273)
Fax: +81-3-5800-6965
E-mail: sada@ua.t.u-tokyo.ac.jp

Accuracy of Areal Interpolation: A Comparison of Alternative Methods
Abstract
This paper discusses the accuracy of spatial data estimated by areal interpolation. A
stochastic model is proposed which represents areal interpolations in diverse geographic
situations. The model is used to examine the relationship between estimation accuracy
and the spatial distribution of estimation error from a theoretical viewpoint. The analysis
shows that the uniformity in error distribution improves the accuracy of areal
interpolation. Four areal interpolation methods are then assessed through numericalexaminations. From this it is found that the accuracy of simple interpolation methods
heavily depends on the appropriateness of their hypothetical distributions, whereas theaccuracy of intelligent methods depends on the fitness of the range of supplementary data
for that of true distribution.
Keywords: areal interpolation; accuracy; error distribution; stochastic model.

1. INTRODUCTION
Spatial data, especially socioeconomic data, are often provided in an aggregatedform though they are originally disaggregated. Census data, for instance, are aggregatedacross census tracts in order to keep confidentiality of subjects. Landuse data aresometimes aggregated across administrative units or square lattice for the reduction ofdata volume. Since there are a wide variety of zonal systems used for data aggregation, it
is a basic function of GIS to integrate multiple spatial databases having incompatible zonal
systems into one database on a common zonal system.
Integration of spatial data requires a data transfer from one zonal system to another.
This process is called areal interpolation, which is inevitably uncertain to some extent
because it involves data estimation based on arbitrary assumptions on the distribution of
spatial objects. Hence, in order to improve estimation accuracy, numerous areal
interpolation methods have been proposed in the literature (Wright, 1936; Markoff andShapiro, 1973; Tobler, 1979; Goodchild and Lam, 1980; Lam, 1983; Flowerdew andGreen, 1991; Rhind, 1991; Goodchild et al., 1993; Burrough and McDonnell, 1998).
Along with the development of new methods, there has arisen a need forcomparison of estimation accuracy between interpolation methods. To meet this demand,
geographers have recently compared the accuracy of areal interpolation methods.
Langford et al. (1991) investigated the accuracy of three areal interpolation methods.
They estimated the population distribution of northern Leicestershire by regression
models where landuse data were used as independent variables. Goodchild et al. (1993)
proposed an interpolation method using ancillary data called "control zones" in which
population density was reasonably assumed to be constant. They employed seveninterpolation methods including their new method to estimate the population of counties in
California and compared the accuracy of estimates.
Most of comparative studies including the above two papers have been performed
on a single geographic situation. Hence it is questionable whether the obtained results
have global applicability (Fisher and Langford 1995). To solve the problem, Fisher,
Langford, and their colleagues have advocated the use of Monte Carlo simulation(Langford et al., 1993; Fisher and Langford, 1995, 1996; Cockings et al., 1997). Monte
Carlo simulation allows us to test interpolation methods in a variety of geographic
circumstances so that the results have wider applicability. Fisher and Langford (1995),
for instance, discussed estimation accuracy of five interpolation methods considering
diverse combinations of source and target zonal systems. They successfully obtained full
distributions of estimation errors which were valuable information to assess theinterpolation methods.
Unfortunately, Monte Carlo simulation is computationally expensive. The

1computational cost is a serious problem especially in its application to areal interpolation,
because it involves the polygon overlay operation which requires rather complicatedalgorithms in GIS. Realization of a number of spatial relationships between source andtarget zones requires considerable calculation. Furthermore, estimation accuracy depends
not only on the relationship between source and target zones but also on the distributionof spatial objects, say, population distribution. Analysis of its effect imposes further
complexity on computation and consequently makes Monte Carlo simulation impractical.
Motivated from the above discussion, we employ a theory-based approach tocompare estimation accuracy among areal interpolation methods. Our approach is basedon a stochastic model that represents diverse geographic situations, and the model
requires far less computational cost than Monte Carlo simulation. In addition to thisadvantage, it allows us to discuss theoretically the relationship between estimation
accuracy and areal interpolation methods.
In the following section we outline the process of areal interpolation and propose amethod for analyzing its accuracy. In Section 3 we discuss the relationship between theaccuracy of areal interpolation and the spatial distribution of estimation errors from atheoretical point of view. In Section 4 we numerically examine the accuracy of four areal
interpolation methods. Finally we summarize the conclusions obtained in this paper withdiscussion.
2. AREAL INTERPOLATION MODEL
Areal interpolation is a process of transferring data from one zonal system toanother. Let us consider, for instance, a zonal system shown in Figure 1a. Point objects
such as individual people or households are distributed over the zones, and the number ofpoint objects is recorded in each zone. A target zone is overlaid on the zones (Figure 1b)
and we want to know the number of points contained in the target zone which is directly
unobservable. In such a case we employ areal interpolation to estimate the number ofpoints from the source data.
Figure 1. Areal interpolation. a) A zonal system, b) an overlaid target zone. The number
in parenthesis indicates the number of points in each zone.
Areal interpolation can be mathematically represented as follows. In this paper we
focus on the estimation of count data in the union of a source zone S and a target zone T.
Let f(x) be the density function of point objects in S. In usual areal interpolation methods
an estimator function of f(x) is implicitly employed. For instance, the areal weighting
method assumes a constant density in S, that is, it assumes f(x) to be a uniform2

distribution. Pycnophylactic interpolation method uses another smooth continuous
surface (Tobler, 1979). Hence we denote an estimator of f(x) as f̂ ( x ) which we callestimator function (Figure 2). Note that the volume-preserving principle (Lam, 1983)
restricts f̂ ( x ) to satisfy

∫

x ∈S

f̂ ( x )dx = ∫x ∈Sf ( x )dx .

(1)

For convenience we denote the volume shown above as V hereafter.
Figure 2. A density function f(x) and its estimator function f̂ ( x ) in one dimensional case.
In the above setting, the number of points in S∩T and its estimator are given byM = ∫ f ( x )1( x )dx
(2)
x ∈S

and
ˆ =
M
∫

x ∈S

fˆ ( x )1( x )dx ,(3)

respectively, where 1(x) is the indicator function of T defined by1 if x ∈ T.
1( x ) = 
0 otherwise(4)

To evaluate the accuracy of areal interpolation, we consider a variety of spatial
relationships between the source and target zones. More explicitly, we assume that thetarget zone T is dropped randomly in such a way that it intersects S, given the density
function f(x), its estimator f̂ ( x ) , and the shape and size of S and T (similar approach was
employed in Sadahiro, 1998a, 1999b). We call this stochastic process areal interpolation
model. The accuracy of interpolation is measured by the mean square error given by2
2
E  M̂ − M  = E  ∫ ∆f ( x )1( x )dx  ,(5)
 x ∈Swhere(6)∆f ( x ) = fˆ ( x ) − f ( x ) .

(

)

{

}

We refer the function ∆f(x) as error distribution function since it represents the spatial
distribution of estimation errors. After a few steps of calculation (see Sadahiro 1998b for
details), we obtain2
1
ˆ −M =E  M(7)m(T ; x − t )∆f ( x )∆f (t )dxdt , m' (T ; S ∩ T ≠ ∅) ∫t ∈S ∫x ∈S()

where m(T; S∩T≠∅) and m(T; l) are the measure of the set of all figures congruent to T
intersecting S and that containing two points separated by a distance l, respectively(Santaló, 1976). These measures are computable by at least numerical calculations(Sadahiro, 1999a), hence it is possible to compare areal interpolation methods in terms of
estimation accuracy using equation (7).
The above equation, however, is somewhat intractable because the computation of3

m(T; l) often involves complicated geometrical calculations. To obtain a more operational
representation of accuracy, we assume that the target zone T is large enough that theincircle of T is larger than the circumcircle of S. If both S and T have rounded shapes, thesize of these zones can be almost equal. This assumption seems realistic in areal
interpolation, because S and T are usually convex in a global scale or T is fairly larger
than S (otherwise areal interpolation may yield considerably imprecise estimates).
If the above assumption holds, we can approximate m(T; l) bym(T ; l ) ≈ 2πAT − 2 LT l ,(8)where AT and LT are the area and perimeter of T, respectively (Stoyan and Stoyan, 1994;
Sadahiro, 1998a). Substitution of equation (8) into equation (7) yields
22 LTˆ −M =−E  M(9)x − t ∆f ( x )∆f (t )dxdt .
m' (T ; S ∩ T ≠ ∅) ∫t ∈S ∫x ∈S()

Since we are interested in comparing various areal interpolation methods, we keep the
shape and size of S and T fixed hereafter. Then it is enough to considerIε = − ∫ ∫ x − t ∆f ( x )∆f (t )dxdt(10)t ∈S x ∈S

for the purpose of assessing the accuracy of areal interpolation. We call this value error
index, and use it as a measure of estimation accuracy in the following sections.
3. ACCURACY OF AREAL INTERPOLATION AND SPATIAL
DISTRIBUTION OF ESTIMATION ERRORSAs seen in equation (10), the accuracy of areal interpolation depends on the error
distribution represented by ∆f(x). To understand this relationship more deeply, weinvestigate how the form of ∆f(x) affects the error index Iε taking some typical examples.
In the following discussion, we focus on a small subregion of S∩T denoted by U
(Figure 3). Keeping ∆f(x) in S∩T\U fixed, we discuss how the form of ∆f(x) in U
affects Iε. The region U is divided into four congruent subregions (not necessarilysquares) U1, U2, U3, and U4, in each of which ∆f(x) has a constant value (Figure 4).
Figure 3 The region U in S∩T.
Figure 4 An example of the error distribution ∆f(x) in U.
3.1 Degree of Estimation Errors
We start with the two cases shown in Figure 5. As shown in the figure, theabsolute value of ∆f(x) in Case 1b is larger than that in Case 1a by d (d>0).
Figure 5 The error distribution function ∆f(x) in U (a>0, b<0).

4

We wish to examine in which case estimates are more accurate. Though the answer
is intuitively obvious, it would be meaningful to discuss this question in a theoreticalway.
Let us denote the error index of the Case i as Iε(i). The index Iε(1a) is then writtenasIε (1a ) = − ∫∫t ∈S x ∈S= −∫x − t ∆f ( x )∆f (t )dxdt∫

t ∈S \ U x ∈S \ U

x − t ∆f ( x )∆f (t )dxdt−2 a( K3 + K4 ) − 2b( K1 + K2 )

,

(11)

− a 2 (k33 + 2 k34 + k44 ) − b 2 (k11 + 2 k12 + k22 )−2 ab(k13 + k14 + k23 + k24 )where
Ki = ∫∫

x − t dxdt

(12)x − t dxdt .

(13)t ∈S \ U x ∈U i

andkij = ∫t ∈U i∫x ∈U j

The index Iε(1b) is also given by a similar equation. Then we haveIε (1b) − Iε (1a ) = 2 dK1 + 2 dK2 − 2 dK3 − 2 dK4+4 d (b − a)(k11 + k12 )−4 d (b − a − d )(k13 + k14 )

.

(14)

≈ 4 d {−( a − b)(k11 + k12 ) + ( a − b + d )(k13 + k14 )}= 4 d {( a − b)(k14 − k11 ) + d (k13 + k14 )} > 0
From equation (14) we obtainIε (1b) > Iε (1a ) .

(15)This confirms that estimation accuracy increases as the absolute value of the error
distribution function ∆f(x) decreases.
3.2 Uniformity in Error Distribution
We then proceed to the second example shown in Figure 6. Though the sum of∆f(x) in U is equal in both cases, its spatial allocation is different: in Case 2a the error is
uniformly distributed whereas in Case 2b the error distribution has spatial variation.
Using this example we investigate whether the uniformity of ∆f(x) improves estimation
accuracy.
Figure 6 The error distribution function ∆f(x) in U.

5

The difference between Iε(2a) and Iε(2b) is given byIε (2 b) − Iε (2a ) = 2 dK1 + 2 dK2 − 2 dK3 − 2 dK4+4 d 2 (k14 − k11 )

.

(16)

≈ 4 d 2 (k14 − k11 ) > 0
From equation (16) we haveIε (2 b) > Iε (2a ) .

(17)

Inequality (17) implies that the spatial variation of estimation error increases the error
index Iε. This leads to a conclusion that the uniformity of ∆f(x) improves the accuracy of
areal interpolation.
3.3 Scale of Uniformity in Error Distribution
The analysis of Case 2 suggests that the uniformity is a desirable property of ∆f(x).
The uniformity, however, can be obtained in a variety of scales, from global to local.
Consequently, we next examine which uniformity in ∆f(x) is more desirable, global(Case 3a) or local (Case 3b).
Figure 7 The error distribution function ∆f(x) in U.
The index Iε(3a) is written asIε (3a ) = − ∫∫t ∈S x ∈S= −∫∫x − t ∆f ( x )∆f (t )dxdt

t ∈S \ U x ∈S \ U

x − t ∆f ( x )∆f (t )dxdt−2 a( K2 + K3 ) − 2b( K1 + K4 )

.

(18)

− a 2 (k11 + 2 k14 + k44 ) − b 2 (k22 + 2 k23 + k33 )−2 ab(k12 + k13 + k24 + k34 )The index Iε(3b) is also given by a similar equation. The difference of the indices is givenbyIε (3b) − Iε (3a ) = 2( a − b)( K2 − K4 )+2(k14 − k12 )( a 2 − 2 ab + b 2 ) .

≈ 2(k14 − k12 )( a − b)If a≠b then equation (19) is positive. Thus we haveIε (3b) > Iε (3a ) .

(19)

2

(20)

Inequality (20) suggests that the global uniformity is more desirable than the local
uniformity. Areal interpolation should be performed so that f̂ ( x ) fits globally for f(x).
From the above analyses, we obtained some desirable properties of areal

6interpolation methods. The first analysis showed that the estimation accuracy increases as
the absolute value of the error distribution function ∆f(x) decreases. This is an intuitively
understandable result because Iε increases in proportion to the square of ∆f(x) as shownin equation (10). In the second analysis we found that the uniformity of ∆f(x) improves
estimation accuracy. This implies that local but serious estimation error may considerably
reduce estimation accuracy. The third analysis showed that the global uniformity is more
important than the local uniformity in terms of estimation accuracy. These results can besummarized as shown in Figure 8.
Figure 8 The error distribution function ∆f(x) and estimation accuracy.
4. NUMERICAL EXAMINATIONSIn the previous section we have investigated the relationship between the accuracyof areal interpolation and spatial distribution of estimation errors. Since we employed atheory-based method in the analysis, the obtained results are rigid and deepen our
understanding of areal interpolation. On the other hand, since no specific interpolation
method was referred in the discussion, one might think that the results are abstract anddifficult to apply to the choice of areal interpolation methods. Hence in this section weanalyze the accuracy of areal interpolation in a more concrete way. To this end, we
numerically evaluate four areal interpolation methods that are widely used in geography
and GIS using the areal interpolation model proposed in Section 2.
We adopt a circle of radius 1.0 as the source zone S. This is because the circle is agood approximation of convex figures in terms of the accuracy of areal interpolation(Sadahiro 1999b). The representative point is located at the center of S whose locational
vector is indicated by z (Figure 9). For the volume V, we try both 1/2π and 2π.
Figure 9 The source zone S and its representative point.
4.1 Areal Interpolation Methods
In numerical examinations we consider four methods of areal interpolation: 1) areal
weighting, 2) point-in-polygon, 3) kernel, and 4) intelligent methods. The first threemethods are often called simple methods in contrast to intelligent methods which use
supplementary data in estimation process. We successively outline these methods with therepresentation of the estimator function f̂ ( x ) (for derivations of f̂ ( x ) , see Sadahiro,
1998b).
(1) Areal weighting method

7

The areal weighting method is one of the most popular simple methods, which isfrequently used when no information is available on the distribution of spatial objects.
This method assumes that spatial objects are uniformly distributed in the source zone
(Figure 10a). Consequently, the estimator function f̂ ( x ) is written asf̂ ( x ) =

∫

x ∈S

f ( x )dx

π

=

V
.
π

(21)

(2) Point-in-polygon method
In GIS the number of spatial objects in a zone is often recorded as an attribute of therepresentative point of the zone. The point-in-polygon method sums up all the counts
allocated to representative points that are included in the target zone (Okabe and Sadahiro,
1997; Sadahiro, 1998a). This implies that the method assumes all the spatial objects to belocated exactly on the representative point (Figure 10b). Hence the estimator functionf̂ ( x ) is given by a delta function which is represented as the equation
 V
ˆf ( x ) =  πr 2 x ∈ C(z, r )(22)
 0 x ∉ C(z, r )in the limit r→0, where C(z, r) is the circle of radius r whose center is located at z.
(3) Kernel method
Kernel method is originally a statistical method used for nonparametric density
estimation (Silverman, 1986; Scott, 1992). However, it is also applicable to areal
interpolation when the locational data of representative points are available (Bracken andMartin, 1989, 1995; Bracken, 1993). Kernel method treats the representative point as a"high information point," that is, it assumes point objects to be clustered around therepresentative point. The probability density function called kernel is placed on the
representative points to create a smooth surface.
There have been proposed various kernel functions in the literature. From thesefunctions we choose the conical kernel for numerical examinations (Figure 10c).
Mathematical representation of f̂ ( x ) is somewhat complicated because of the volumepreserving principle. Suppose that the slope of the cone λ is given. The estimatorfunction f̂ ( x ) is then represented as follows:1
if πλ < V then
3
1
2
(23)
f̂ ( x ) =  V + πλ  − λ x − z ,π
3 otherwise,8 3λ2V3f̂ ( x ) =  π − λ x − z0
For the value of λ we try 1.0, 3.0 and 5.0.

3Vπλ .
otherwise

if x − z ≤ 3(24)(4) Intelligent method
As GIS and remotely-sensed satellite images become widely available, new areal
interpolation methods called intelligent methods have been developed (Fisher andLangford, 1995, 1996; Langford et al., 1991; Goodchild et al., 1993). They use
supplementary data such as satellite images or landuse data in areal interpolation to
improve estimation accuracy.
Intelligent methods can be described mathematically by use of a supplementaryfunction ϕ(x) which represents additional data. Let us consider, for instance, thedasymetric method which is often used for population estimation. Using landuse data thismethod divides a source zone into subregions in each of which population density is
assumed to be constant. Typically, a source zone is divided into two subregions, that is,
residential and non-residential areas. In this case the function ϕ(x) is written as a binary
function
1 if x ∈ R
ϕ (x) = ,(25)
0 otherwise
where R indicates the residential area. The estimator function f̂ ( x ) is then given byϕ (x)f̂ ( x ) =V.
(26)∫ ϕ (x)dxx ∈S

Since population density is highly correlated with landuse category, equation (25) may be
represented as
1 if f ( x ) ≥ α,(27)
ϕ (x) = 
0 otherwise
where α is the lower limit of the population density of residential areas.
If more than two landuse categories are distinguishable, ϕ(x) is given by a stepwisefunctionϕ ( x ) = ki if x ∈ Li
(28)
(Figure 10d), where Li indicates the type i landuse area. Substituting equation (28) into
equation (26), we obtain the estimator function f̂ ( x ) .
Linear regression models can be regarded as the case where the ϕ(x) is a linear
function of f(x) having a probabilistic error term ε, that is,ϕ ( x ) = α + βf ( x ) + ε .
(29)Non-linear regression models can also be represented in a similar way.

9In numerical examinations, we try three forms of the supplementary function ϕ(x)as follows.
1 if f ( x ) ≥ 1(30)
ϕ1 ( x ) = 
Method I1:
0 otherwise
if f ( x ) ≥ 1
 1
2 / 3 if 1 > f ( x ) ≥ 2 / 3Method I2:
ϕ 2 (x) = 
(31)
1
3
2
3
>
x
≥
1
3
/
if
/
f
/
(
)

 0
if 1 / 3 > f ( x )if f ( x ) ≥ 3
3
2 if 3 > f ( x ) ≥ 2

ϕ 3 (x) = 
Method I3:
(32)
1
2
>
x
≥
1
if
f
(
)

0
if 1 > f ( x )If ϕ(x)=0 for any x, we redefine the function as
ϕ (x) = k .

(33)

In this case the intelligent method is equivalent to the areal weighting method.
Figure 10 Estimator functions of areal interpolation methods. a) Areal weighting, b)
point-in-polygon, c) kernel, and d) intelligent methods.
4.2 Results
In order to assess the accuracy of the above interpolation methods, we employ fourtypes of point distributions as f(x): 1) cylinder, 2) conical, 3) annular, and 4) linear
distributions (Figure 11). The results of numerical examinations are successively
described in the following.
Figure 11 Forms of the function f(x). a) Cylinder, b) conical, c) annular, and d) linear
distributions.
(1) Cylinder distribution
Equation (10) and the results obtained in Section 3 indicate that areal interpolation
yields better estimates if the estimator function is close to the true distribution.
Consequently, for instance, the point-in-polygon and kernel methods are expected to give
better results when point objects are clustered around the representative point. To confirm
this we examine how the degree of concentration affects the accuracy of estimates using
the cylinder distribution, a typical concentrated distribution of point objects. The radius of
the cylinder indicated by r changes from 0 to 1 as shown in Figure 12. The height of thecylinder is denoted by h.

10Figure 12 Cylinder distributions.
The relationship between the radius r and estimation accuracy is depicted in Figure
13. The results are quite consistent with our intuition. As the distribution becomes
concentrated estimation accuracy of the point-in-polygon method linearly increases
whereas that of the areal weighting method decreases. The kernel methods are in the
middle of these methods. The accuracy of the intelligent methods is noteworthy: the error
index Iε is nearly always zero. This is because the stepwise functions assumed in theintelligent methods fit locally-uniform distributions.
Figure 13 Estimation accuracy on the cylinder distribution. r: the radius of the cylinder. h:
the height of the cylinder. Ki : the kernel method of λ=i. Ii : the intelligent method of thefunction ϕi (x). a) V=1/2π, b) V=2π.
(2) Conical distribution
We next discuss the conical distribution as another example of centralized
distributions. The slope of the cone denoted by m increases gradually from 0 to 100(Figure 14). The height of the cone is indicated by h.
Figure 14 Conical distributions.
The results are shown in Figure 15. Concerning the areal weighting and point-inpolygon methods the results are quite understandable as in the case of the cylinder
distribution. On the other hand, it is counterintuitive that the intelligent methods do notalways give the best result among the four interpolation methods. For instance, let us takea look at Figure 15a-2. When the slope m is between 1 and 2, the intelligent methods I1
and I3 are far less accurate than the areal weighting and kernel methods. The intelligent
methods do not suit the gently-sloping distribution, though they work successfully for the
distribution of a steep slope (see Figures 15a-1 and 15b-1). This can be partly understoodby looking carefully at Figure 15b-2. When the cone height h is around 5.0, theintelligent method I3 whose boundary values are 1.0, 2.0, and 3.0 (recall equation (32))
gives better estimates than I1 and I2 having boundaries of smaller values. For h=2.1,
however, I3 is the least accurate among the intelligent methods. From this we can say that
estimation accuracy of the intelligent methods depends on the fitness of the range of the
stepwise function ϕ(x) for that of the density function f(x). They yield poor estimates if
the boundaries of ϕ(x) do not agree with the range of f(x).

11Another counterintuitive result is also found in Figure 15b. The intelligent method
I2 is less accurate than I1 for h>4.0, though I2 has a finer stepwise function than I1. Thisis also caused by the lack of agreement between the ranges of ϕ(x) and f(x). The finenessof the stepwise function sometimes makes areal interpolation less accurate, if the rangesof supplementary data and true distribution do not agree.
One might think it strange that the intelligent methods yield the best estimates forvery large h, say, h=20.0. In such a case, however, f(x) is equal to zero in almost
everywhere so that the boundary setting of ϕ(x) scarcely affects the estimation accuracy.
Figure 15 Estimation accuracy on the conical distribution. m: the slope of the cone. h: theheight of the cone. Ki : the kernel method of λ=i. Ii : the intelligent method of the function
ϕi (x). a-1) V=1/2π, a-2) V=1/2π (magnified), b-1) V=2π, b-2) V=2π (magnified).
(3) Annular distribution
Contrary to the previous two distributions we next consider a disconcentrated
distribution, that is, the annular distribution. Using this distribution we analyze therelationship between the degree of discentralization and the estimation accuracy of areal
interpolation. The radius of the hole indicated by r changes from 0 to 1 as shown in
Figure 16. The height of the annulus is denoted by h.
Figure 16 Annular distributions.
The results are shown in Figure 17. It depicts that the simple methods fail to giveaccurate estimates when point objects are distributed along the edge of the source zone.
This is obviously because the annular distribution is quite different from the distribution
assumed in the simple methods. Only the intelligent methods yield satisfactory results forthe annular distribution. The results warn us that it is quite dangerous to assume blindly
an arbitrary distribution in areal interpolation.
Figure 17 Estimation accuracy on the annular distribution. r: the radius of the hole. h: theheight of the annulus. Ki : the kernel method of λ=i. Ii : the intelligent method of thefunction ϕi (x). a) V=1/2π, b) V=2π.
(4) Linear distribution
We finally examine the linear distribution. This type of monotonous distribution,
though not necessarily linear, is often found in global population distribution. The slope
indicated by m changes from 0 to V/π as shown in Figure 18.

12Figure 18 Linear distributions.
The results shown in Figure 19 are similar to those in Figure 17 in that the simple
methods are generally inaccurate. Estimation accuracy of the simple methods
monotonously reduces with an increase of m. On the other hand, the intelligent methods
yield better estimates for any m. Among the intelligent methods, I1 and I2 are better thanI3 for h<2.6 (Figure 19b), whereas I3 is better for larger h. This tendency was also foundin the conical distribution case (Figure 15b), which is caused by the disagreement
between the ranges of ϕ(x) and f(x) as discussed earlier.
Figure 19 Estimation accuracy on the linear distribution. m: the slope of the linear
distribution. h: the height of the distribution. Ki : the kernel method of λ=i. Ii : theintelligent method of the function ϕi (x). a) V=1/2π, b) V=2π.
5. CONCLUDING DISCUSSION
In this paper we have analyzed the accuracy of areal interpolation from boththeoretical and empirical viewpoints. To keep generality in the analysis, we firstdeveloped a stochastic process which we call areal interpolation model that represents
areal interpolation in diverse geographical situations. We then proposed a measure ofestimation accuracy called error index on the basis of the model.
The results obtained from theoretical considerations in Section 3 can be summarized
as follows: 1) estimation accuracy improves as the error at individual locations decreases;2) estimation accuracy improves as the error distribution becomes smooth and uniform; 3)
the global uniformity in the error distribution is more important than the local uniformity
in terms of estimation accuracy. Though the first conclusion is not surprising, the second
and third results are useful for the evaluation of areal interpolation methods. The spatial
uniformity in the error distribution, or the global fitness of estimator function, is a
desirable property of areal interpolation methods.
From numerical examinations we obtained the following results: 1) intelligent
methods are generally more accurate than simple methods; 2) estimation accuracy ofsimple methods heavily depends on the appropriateness of their hypothetical distributions;3) estimation accuracy of intelligent methods depends on the fitness of the range ofsupplementary data for that of true distribution; 4) if the ranges do not agree, intelligent
methods may yield poorer estimates than simple methods.
The above results are quite suggestive for the choice of areal interpolation methods.
However, they are clearly limited in several aspects. We finally discuss the limitations of13this paper for further research.
First, the effects of the shape and size of the target and source zones are remained
to be analyzed. We have fixed these factors throughout the analysis in order to focus on
how the relationship between two spatial distributions, one is actual and the other isassumed in areal interpolation, affects estimation accuracy. Though the influence of thisrelationship is quite significant, geometrical properties of spatial units are also important.
Concerning the areal weighting method, their effects on estimation accuracy werenumerically investigated (Sadahiro, 1999b). However, the study was limited to one
interpolation method, so accuracy of other methods should be analyzed in relation togeometrical properties of zones.
Second, we have discussed the accuracy of areal interpolation only for a pair ofsource and target zones. This allows us to examine the relationship between estimation
accuracy and error distribution from a theoretical point of view. However, there is a clear
demand for theoretical analysis in more realistic situations - a combination of multiple
source and target zones. Though it is quite difficult to discuss estimation accuracy
keeping generality in complicated geographical circumstances, the difficulties have to be
overcome in future research.
Third, there still remain areal interpolation methods not examined in this paper. For
instance, the Poisson regression method (Flowerdew and Green, 1989) and thepycnophylactic method (Tobler, 1979) should also be tested. Since they involve a
convergence calculation in estimation process, the first step to assess their accuracy is toanalyze the properties of the estimator function yielded from convergence calculation.
Once the estimator function is given in an analytical form, their accuracy can be evaluatedby the method proposed in this paper.

14

LITERATURE CITED
Bracken, I. and D. Martin (1989). "The Generation of Spatial Population Distributions
from Census Centroid Data." Environment and Planning A 21, 537-543.
Bracken, I. (1993). "An Extensive Surface Model Database for Population-RelatedInformation." Environment and Planning B 20, 13-27.
Bracken, I. and D. Martin (1995). "Linkage of the 1981 and 1991 UK Censuses Using
Surface Modelling Concepts." Environment and Planning A 27, 379-390.
Burrough, P. A. and R. A. McDonnell (1998). Principles of Geographical Information
Systems. New York: Oxford University Press.
Cockings, S., P. F. Fisher, and M. Langford (1997). "Parameterization and
Visualization of the Errors in Areal Interpolation." Geographical Analysis 29, 314328.
Fisher, P. F. and M. Langford (1995). "Modelling the Errors in Areal Interpolation
between Zonal Systems by Monte Carlo Simulation." Environment and Planning A
27, 211-224.
Fisher, P. F. and M. Langford (1996). "Modeling Sensitivity to Accuracy in Classified
Imagery: A Study of Areal Interpolation by Dasymetric Mapping." Professional
Geographer 48, 299-309.
Flowerdew, R. and M. Green (1989). "Statistical Methods for Inference between
Incompatible Zonal Systems." In Accuracy of Spatial Databases, edited by M.
Goodchild and S. Gopal, pp. 239-247, London: Taylor and Francis.
Flowerdew, R. and M. Green (1991). "Data Integration: Statistical Methods forTransferring Data between Zonal Systems." In Handling Geographical Information:
Methodology and Potential Applications, edited by I. Masser and M. Blakemore, pp.
38-54, New York: Longman.
Goodchild, M. F. and N. N-S. Lam (1980). "Areal Interpolation: a Variant of theTraditional Spatial Problem." Geo-processing 1, 297-312.
Goodchild, M. F., L. Anselin and U. Deichmann (1993). "A Framework for the Areal
Interpolation of Socioeconomic Data." Environment and Planning A 25, 383-397.
Lam, N. N-S. (1983). "Spatial Interpolation Methods: a Review." American
Cartographer 10, 129-149.
Langford, M., D. J. Maguire, and D. J. Unwin (1991). "The Areal Interpolation
Problem: Estimating Population Using Remote Sensing in a GIS Framework." In
Handling Geographical Information: Methodology and Potential Applications,edited by I. Masser and M. Blakemore, pp. 55-77, New York: Longman.
Langford, M., P. Fisher, and D. Troughear (1993). "Comparative Accuracy
Measurements of the Cross Areal Interpolation of Population." Proceedings of EGIS

15

'93, 663-674.
Markoff, J. and G. Shapiro (1973). "The Linkage of Data Describing Overlapping
Geographical Units." Historical Methods Newsletter 7, 34-46.
Okabe, A. and Y. Sadahiro (1997). "Variation in Count Data Transferred from a Set of
Irregular Zones to a Set of Regular Zones through the Point-in-polygon Method."International Journal of Geographical Information Science 11, 93-106.
Rhind, D. W. (1991). "Counting the People: the Role of GIS." In Geographical
Information Systems, Volume 2: Principles and Applications, edited by D. J.
Maguire, M. F. Goodchild, and D. W. Rhind, pp. 127-137, New York: Longman.
Sadahiro, Y. (1998a). "Accuracy Count Data Estimated by the Point-in-Polygon
Method." Discussion Paper Series 77E, Department of Urban Engineering,
University of Tokyo.
Sadahiro, Y. (1998b). "Accuracy of Areal Interpolation for Surface." Discussion Paper
Series 79, Department of Urban Engineering, University of Tokyo (in Japanese).
Sadahiro, Y. (1999a). "Statistical Methods for Analyzing the Distribution of Spatial
Objects in Relation to a Surface." Journal of Geographical Systems, to appear.
Sadahiro, Y. (1999b). "Accuracy of Count Data Transferred through the Areal Weighting
Interpolation Method." International Journal of Geographical Information Science,
to appear.
Santaló, L. A. (1976). Integral Geometry and Geometric Probability. London: AddisonWesley.
Scott, D. W. (1992). Multivariate Density Estimation. New York: John Wiley.
Silverman, B. W. (1986). Density Estimation for Statistics and Data Analysis. London:
Chapman & Hall.
Stoyan, D. and H. Stoyan (1994). Fractals, Random Shapes and Point Fields. New
York: John Wiley.
Tobler, W. R. (1979). "Smooth Pycnophylactic Interpolation for Geographical Regions."
Journal of the American Statistical Association 74, 519-530.
Wright, J. K. (1936). "A Method of Mapping Densities of Population with Cape Cod as
an Example." Geographical Review 26, 103-110.

16

zone 1(17)
zone 4(25)

zone 2(10)zone 5(17)target zone(?)

zone 3(31)
zone 6
(9)

(a)

(b)

Figure 1

density function f(x)estimator function f(x)Figure 2

T

US

Figure 3∆f(x)

∆f(x)

U3
0

U1

U4
U2

0

Figure 4∆f(x)

∆f(x)

a+d
a

0

0

b
b-d

Case 1aCase 1b

Figure 5∆f(x)

∆f(x)

a+d
a
a-d
0

0

Case 2aCase 2b

Figure 6∆f(x)

∆f(x)

a

a

b
0

b
0

Case 3aCase 3b

Figure 7

high

low
estimation accuracy

Figure 8z

1.0

Figure 9

(a)

(b)

(c)

Figure 10

(d)

1.0
(a)

(b)

(c)

Figure 11

(d)

h

r

1.0

Figure 12

1.2

AW: Areal weighting method
PIP: Point-in-polygon method
Ki: Kernel methods
Ii: Intelligent methods

1.0
PIP
0.8

AW
Iε 0.6

K1
K3

0.4

K5
0.2
I1, I3
I2

0.0
1.0
0.5

0.8
0.6

0.8

0.6
1.0

1.4

0.4
2.0

Figure 13a

3.1

0.0 r

0.2
5.6

12.5

50.0

∞ h

18
16

AW: Areal weighting method
PIP: Point-in-polygon method
Ki: Kernel methods
Ii: Intelligent methods

14
PIP
12

AW

10

K1

Iε

K3

8

K5
6
4
2
I1, I2, I3

0
1.0
2.0

0.8
2.5

3.1

0.6
4.1

5.6

0.4
8.0

12.5

Figure 13b

0.0 r

0.2
22.2

50.0

200.0

∞ h

h

m

r

1.0

Figure 14

1.2

AW: Areal weighting method
PIP: Point-in-polygon method
Ki: Kernel methods
Ii: Intelligent methods

1.0

0.8
AW
Iε 0.6
PIP

K1

0.4
K3
0.2
K5
I1 I2 I3

0.0
0
0.5

20
5.3

8.4

40
11.1

13.4

60
15.5

17.5

Figure 15a-1

100 m

80
19.4

21.3

23.0

24.7 h0.25
I1, I3
AW: Areal weighting method
PIP: Point-in-polygon method
Ki: Kernel methods
Ii: Intelligent methods

0.20

AW
0.15
K5
Iε
0.10
K3

0.05
K1

I3

I1

I2

0.00
0
0.5

1
0.8

1.2

2
1.5

1.8

3
2.1

2.4

Figure 15a-2

5 m

4
2.6

2.9

3.1

3.3 h

18
16

AW: Areal weighting method
PIP: Point-in-polygon method
Ki: Kernel methods
Ii: Intelligent methods

14
12
10
Iε

PIP
8

AW

6

K1

4

K3
K5

2
I1 I2 I3
0
0
2.0

20
8.4

13.4

40
17.5

21.3

60
24.7

27.8

Figure 15b-1

100 m

80
30.9

33.7

36.5

39.1 h0.7

AW: Areal weighting method
PIP: Point-in-polygon method
Ki: Kernel methods
Ii: Intelligent methods

0.6

0.5
AW
0.4
K5

Iε
0.3

I2
0.2
K3
I1
0.1
K1
0.0

I3

0
2.0

1
2.3

2.7

2
3.0

3.3

3
3.7

4.0

Figure 15b-2

5 m

4
4.3

4.7

5.0

5.3 h

h

r

1.0

Figure 16

1.8
1.6
1.4
1.2
AW: Areal weighting method
PIP: Point-in-polygon method
Ki: Kernel methods
Ii: Intelligent methods

PIP
1.0
Iε
0.8
0.6

K5
K3

0.4

K1

0.2
0.0
0.0
0.5

0.2
0.51

0.52

0.4
0.55

0.59

0.6
0.67

0.78

Figure 17a

I2
0.98

I1, I3

AW
1.0 r

0.8
1.39

2.63

∞ h

30
PIP
25

20
AW: Areal weighting method
PIP: Point-in-polygon method
Ki: Kernel methods
Ii: Intelligent methods

Iε 15

10

5
K5
K3

K1

AW

I 1, I 2, I 3

0
0.0
2.00

0.2
2.02

2.08

0.4
2.20

2.38

0.6
2.67

3.13

Figure 17b

1.0 r

0.8
3.92

5.56

10.53

∞ h

m

h

1.0

Figure 18

1.2
PIP
1.0
AW: Areal weighting method
PIP: Point-in-polygon method
Ki: Kernel methods
Ii: Intelligent methods

0.8

Iε 0.6

0.4
K5
K3
0.2
K1

AW
I2

0.0
0.0
0.5

0.1
0.55

0.6

0.2
0.65

0.7

0.3
0.75

Figure 19a

0.8

0.85

0.4

I1, I3

0.5 m

0.9

0.95

1.0 h

3.5

AW: Areal weighting method
PIP: Point-in-polygon method
Ki: Kernel methods
Ii: Intelligent methods

3.0

2.5

2.0
K5

Iε

K3

K1

AW

1.5
I2

1.0

I1

0.5

I3
0.0

0
2.0

0.4
2.2

2.4

0.8
2.6

2.8

1.2
3.0

Figure 19b

3.2

2.0 m

1.6
3.4

3.6

3.8

4.0 h