CSIS Discussion Paper No. 6
Accuracy of Count Data Estimated by
the Point-in-Polygon Method
Yukio Sadahiro
JANUARY, 1999

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

Accuracy of Count Data Estimated by the Point-in-Polygon Method
Abstract
This paper analyzed the accuracy of count data estimated by the point-in-polygon
method. A point-in-polygon interpolation model was proposed which was based on astochastic distribution of points and the target zone, in order to represent a variety ofsituations. The accuracy of estimates was numerically investigated in relation to the sizeof the target zone and the distribution of points, and the optimal location of representative
points was discussed. The major findings obtained in this paper were as follows: 1)
though the relative accuracy of estimates generally increases monotonously with the sizeof the target zone, the monotoneity is often disturbed by the periodicity in the spatialconfiguration of source zones and the point distribution; 2) the point-in-polygon and theareal weighting interpolation methods have the same accuracy of estimates when points
are concentrated in less than 12-15% area around the representative point in source zones;
3) the point-in-polygon method is not so robust against the locational gap between points
and the representative point; 4) the optimal location of representative points is given by
the spatial median of points.

11. INTRODUCTION
Socioeconomic data based on the point location are often provided in an aggregatedform. Census data, for instance, are aggregated across census tracts to open to the public,
though they are originally collected as point data. Other data are aggregated on a variety ofzonal systems such as regular lattices, municipal districts, school districts, and postal
zones. It often happens that a zonal system is incompatible with the region in which ananalysis is to be performed. In such a case it is necessary to estimate the count data in the
study region. This type of data estimation is called areal interpolation, or to be exact, datatransfer from source zones to a target zone. Since areal interpolation is frequently requiredin geographical analysis, a variety of mathematical methods have been proposed in
geography and other related fields (Wright 1936; Markoff and Shapiro 1973; Tobler
1979; Goodchild and Lam 1980; Lam 1983; Rhind 1991).
Areal interpolation methods can be classified into two groups: intelligent methods
which use supplementary data, say, remotely sensed data, and simple methods which do
not use such data. An advantage of intelligent methods is accuracy of estimates. Using
suitable supplementary data, they usually provide more accurate estimates of spatial data
than simple methods (Flowerdew 1988; Langford et al. 1991; Fisher and Langford
1995). We should note, however, that intelligent methods are not always applicable.
Remotely sensed data are often unavailable or expensive, and the computational cost is
problematic if a great quantity of data has to be processed. Thus simple methods such as
the point-in-polygon method (Burrough and McDonnell 1998), areal weighting
interpolation method (Markoff and Shapiro 1973), and pycnophylactic interpolation
method (Tobler 1979), are still widely used in geography and GIS (Okabe and Sadahiro
1997; Sadahiro 1999b).
Among simple methods, the point-in-polygon method is distinguished for its
processing speed (Flowerdew and Green 1991; Okabe and Sadahiro 1997). This method
transfers the count data of a source zone to a target zone if the representative point of asource zone is included in the target zone (see section 2 for detailed description). The
point-in-polygon method is based on the point-location algorithm (Preparata and Shamos
1985), which runs very fast (linear time). Consequently, the point-in-polygon method issuitable for massive spatial data and in fact it is used for processing 1.65 million census
tracts data by the Japanese Bureau of Statistics (Sinfonica 1994). The importance ofprocessing speed cannot be underestimated in geographical analysis, because in recent
years a number of huge spatial databases have become available to geographers. On theother hand, the point-in-polygon method is said to yield less accurate estimates than othermethods. The method implicitly assumes that all the points are located on representative
points, and gives the true count in a target zone if the assumption holds. However, this2

assumption is apparently too strong. Points are usually dispersed in spatial zones, whichcauses inaccuracy in estimation of count data. Estimation error is relatively small andoften negligible if a target zone is even larger than source zones. Otherwise, estimation
error is crucial to the quality of spatial analysis and thus care should be taken in handlingthe estimated data.
There are two alternatives to improve estimation accuracy when supplementary data
are not available: to use source data consisting of smaller zones, or to adopt anothersimple interpolation method, say, the areal weighting interpolation method. This choice
requires us to understand the nature of estimation accuracy, that is, how accurately areal
interpolation methods estimate count data. Concerning the point-in-polygon method,
unfortunately, there are few systematic studies analyzing the accuracy of estimates.
Though it is known that accuracy depends on various factors such as the size of sourcezones and the location of representative points, their effects have not been quantitatively
evaluated in the literature.
This paper analyzes the accuracy of count data estimated by the point-in-polygon
method in order to investigate how it is determined by various factors. We hope that the
study helps geographers to choose source data and interpolation methods in theiranalysis. In section 2, we briefly describe the point-in-polygon method. We also outline
the areal weighting interpolation method since it is discussed later in comparison with thepoint-in-polygon method. In section 3, we propose a stochastic model representing thepoint-in-polygon interpolation. Using the model, we numerically examine the accuracy of
estimates in section 4 in the case where source zones are a square lattice. In section 5, theoptimal location of representative points is discussed. We finally summarize theconclusions in Section 6.
2. POINT-IN-POLYGON AND AREAL WEIGHTING INTERPOLATIONMETHODS
Suppose a zonal system S, a region Z0 of area A0 consisting of K zones Z1, Z2, ...,ZK (Figure 1a). The region Z0 may represent a ward whereas zone Zi is its lower levelspatial unit, say, a census tract. We call the zones Z 1 , Z 2 , ..., Z K source zones anddenote the area of Zi as Ai (i=1, ..., K). Spatial objects which can be regarded as points,
say, individual people or households, are distributed in Z0 (Figure 1b). We refer to these
spatial objects as points, and denote the location of point j as yj . The number of points iscounted in each source zone, and the data are allocated to the representative points. The
number of points and the location of the representative point in zone Zi are denoted by ni
and pi , respectively (Figure 1c). The locational data of points are then eliminated topreserve the confidentiality of the subjects.

3FIG. 1. A zonal system and count data. a) Source zones, b) distributed points, c)
representative points and the number of points.
We next consider a target zone T of area B in which we want to know the numberof points. Since the location of individual points is not known, we have to estimate thecount using an areal interpolation method. In the point-in-polygon method, the estimate ofcount is given by summing up the ni s whose representative points are included in thetarget zone T (Figure 2a). In the areal weighting interpolation method, the count in eachsource zone is divided according to area and the estimate is given by summing up theassigned values (Figure 2b).
FIG. 2. Areal interpolation methods. a) The point-in-polygon method, b) the areal
weighting interpolation method. The gray-shaded area indicates the target zone. The
numbers in parentheses indicate the count data ni s .
The two methods are mathematically represented as follows. Let us denote thelocation of T using a binary function defined by1 if x ∈ T
C( x ) = .
(1)
0 otherwise
The true number of points in T is then given byM = ∑C yj .

( )

(2)

The number of points in zone Zi is written asni = ∑ Ui y j ,

(3)

j

( )

j

where Ui (x) is defined by1 if x ∈ Zi
Ui ( x ) = .
(4)
0 otherwise
The point-in-polygon method estimates the number of points in T∩Zj as ni if and only if
T contains pi . Otherwise it assigns zero to T∩Zj . Consequently, the estimate of M givenby the point-in-polygon method is written asM̂ = ∑ C(pi )nii( )= ∑ C(pi )∑ Ui y j
i

j

.

(5)

On the other hand, the areal weighting interpolation method estimates the count in T to be

4

∫
ˆ' =
M
∑ x ∈Zii

C( x )dx
Ai

ni .

(6)

3. POINT-IN-POLYGON INTERPOLATION MODEL
Before developing a model, we briefly discuss existing approaches to estimation
accuracy of areal interpolations. The accuracy of areal interpolations has been studied invarious ways, and the approaches taken in the literature can be classified into two
categories: empirical-based studies and Monte Carlo simulations. The former usuallyconsiders a particular geographic situation, that is, actual zonal systems and a point
distribution, in order to evaluate the reliability of estimates and compare the accuracy
between interpolation methods (Flowerdew 1988; Flowerdew and Green 1991; Langford
et al. 1991; Goodchild et al. 1993). Though this approach is advantageous in the reality,
it is an open question whether the obtained results are generally applicable (Fisher andLangford 1995). The Monte Carlo method, on the other hand, enables us to consider awide variety of situations and thus to obtain general results. Fisher and Langford (1995)
and Cockings et al. (1997) adopted this approach to discuss several interpolation methods
in terms of estimation accuracy, and successfully obtained full distributions of estimation
errors. Monte Carlo simulations, however, are computationally expensive, because areal
interpolations involve the polygon overlay operation which requires rather complicatedalgorithms in GIS. The computational cost is problematic especially when a number ofspatial relationships between source and target zones are to be realized.
Considering the discussion above, we follow the approach taken by Sadahiro(1999b), that is, evaluation of estimation accuracy on the basis of a stochastic model
representing a class of situations. This assures us not only the generality in analysis but
also far less computational cost than Monte Carlo simulations.
In the model, a set of source zones with representative points are given whereas
points and the target zone are stochastically located on the source zones. To be explicit,
we consider a situation that N points are independently and identically distributed in theregion Z 0 according to a probability density function f(x), and the target zone T is
dropped in such a way that it intersects Z0. For convenience of computation, we assume
that Z0 has such a shape that can cover a plane by its lattice, and that Z0 is surrounded byits copies having the same zonal system and point distribution as those of Z0 (Figure 3).
This assumption is called the periodic continuation which is often used in spatial statistics(Ripley 1981; Stoyan and Stoyan 1995; Sadahiro 1999b). If T does not completely lie in
Z0, we replace the portion of T outside Z0 by its corresponding figure (Figure 4). The
location of T follows the probabilistic distribution such that all possible shapes and5positions of T appear randomly.
FIG. 3. Periodic continuation assumption.
FIG. 4. Transformation of T.
One might think that the periodic continuation assumption is too strong andunrealistic. The assumption, however, is not essential for the analysis since it isintroduced mainly in order to avoid computational difficulties arising from boundary
effects. Actually the assumption is not necessary if the region Z0 is larger enough than thetarget zone T. There are also several ways to solve boundary problems, and they usually
give almost the same results (Sadahiro 1999a).
In the above setting the accuracy of count data is discussed. The source zonal
system with representative points, the number of points, and the shape and size of thetarget zone are given whereas the location of points and the target zone is probabilisticallydetermined. This implies that we are considering a class of situations that share the given
conditions.
From equations (2) and (5) we have the estimation error of M:
ε = M − Mˆ.
(7)= ∑ C y j − ∑ C(pi )∑ Ui y j

( )

j

i

( )

j

To evaluate the error, we use the mean square error (MSE) of ε defined by
MSE[ S] = E ε 2 .

[ ]

(8)

After several steps of calculation (see Appendix 1 for details) we obtain
MSE[ S] = N ( N − 1)∫t ∈Z 0∫x ∈Z 0f ( x ) f (t ) Pr[ x ∪ t ∈ T ]dxdt

 B
+2 N   A0 
−2 N ( N − 1)∑ ∫x ∈Z ii−2 N ∑ ∫i

x ∈Z if ( x )dx ∫x ∈Z 0f ( x ) Pr[pi ∪ x ∈ T ]dx

f ( x ) Pr[pi ∪ x ∈ T ]dx+ N ( N − 1)∑ ∑ Pr[pi ∪ pi' ∈ T ]∫i

x ∈Z ii'f ( x )dx ∫x ∈Z i '.

(9)f ( x )dx

In the above equation the probability Pr[x∪t∈T] is not explicitly given. It can becomputed as follows. Let m(T; l) be the measure of the set of all figures congruent to T
containing two points separated by a distance l (Santaló 1976). Then Pr[x∪t∈T] is givenby6m( T ; x − u ),2πA0u ∈Ω ( t )Pr[ x ∪ t ∈ T ] = ∑(10)where Ω(t) is a set of points corresponding to t in the surrounding copies of Z0 (Figure5). The measure m(T; l) is represented by an explicit form if T has a simple shape(Santaló 1976; Sadahiro 1999a, 1999b). For instance, if T is a circle of radius r,l4πr 2 arccos  − πl 4r 2 − l 2 (l ≤ 2r ),m( T ; l ) = .
(11) 2r 

0(l > 2r ).
For complicated shapes of T, we can use an equation
B2m( T ; l ) =gT (l ) ,(12)lwhere g T (l) is the probability density function of the distance between two points
randomly distributed in T. Since the function gT(l) is numerically computable, equation(9) can be evaluated for any T.
FIG. 5. A point set Ω(t). The gray-shaded area indicates the region Z0.
4. NUMERICAL EXAMINATIONSHaving obtained a computable representation of MSE, we are now ready to analyze
numerically the accuracy of estimates using the point-in-polygon interpolation model. We
wish to investigate how various factors, especially the size of the target zone and the
distribution of points, affect the accuracy of estimated count data.
As seen in the previous section, the model is applicable for evaluating estimation
accuracy in a variety of situations. However, due to limitations of space, we focus on a
few typical cases where source zones are a square lattice. Though a more realistic zonal
system such as census tracts is also available, we choose a square lattice because of thefollowing reasons. First, congruity of source zones permits us to ignore the effect ofdiversity among source zones on estimation accuracy, and consequently enables us toconcentrate on other specific factors. Second, computational cost is reduced if sourcezones are regularly arranged. This brings high tractability to the numerical examination.
Third, a square lattice is a good approximation of a zonal system whose zones are similarin shape and size. Existing studies such as Okabe and Sadahiro (1997) and Sadahiro(1999b) suggest that a slight difference in shape of source zones does not greatly
influence the accuracy of estimates.
As mentioned above, we focus on how the size of the target zone and the
distribution of points affect the accuracy of estimated count data. To this end, we supposea square lattice whose every cell has its representative point on the centroid of the cell and7has the same form of f(x) (Figure 6). We then consider the limit where the region Z0
expands infinitely, keeping the size of cells and the density of points at A and µ=N/A0,
respectively.
FIG. 6. Representative points and f(x) on a square lattice.
To make the MSE comparable among different sizes of T, we divide it by B2 forstandardization. Hence we have
 1
A f x A f t m T ; x − t )dxdt 
2 ∫t ∈U ∫x ∈Z { 0 ( )}{ 0 ( )} (0
 2πAB
MSE[ S] 2
 1lim= lim − 2 ∫ { A0 f ( x )}m(T ; p − x )dx
2
µ
x
∈
Z
Z 0 →∞
Z
→∞
0
0
B
πB
 , (13)
A2
+ 2πAB2 ∑ m(T ; p − pi )
i


2
1
+  −
A f x m T ; p − x )dx  µ
2 ∫x ∈U { 0 ( )} (
π
B
AB


where U and p indicate a unit cell and a representative point, respectively.
Equation (13) indicates that the mean square error of estimates is represented by alinear combination of µ and µ2. Hence we rewrite the equation as
MSE[ S]lim= Q2 µ 2 + Q1µ ,(14)2Z 0 →∞
Band employ Q1 and Q2 instead of the MSE to evaluate the accuracy of estimates, in order
to discuss the effect of µ separately from other factors. Note that both Q 1 and Q 2 are
independent of µ since they are determined by U, T, p, and f(x).
Given a certain value of µ, the MSE is determined by Q 1 and Q 2. If points aresparsely distributed so that µ has a small value, both Q 1 and Q 2 are influential on
estimation accuracy. However, if points are densely distributed in Z0, we may neglect Q1
in evaluating the MSE since it is mainly governed by Q 2, the leading coefficient inequation (14).
In the following we discuss the point-in-polygon method in comparison with theareal weighting interpolation method because it is also widely used in geography. The
standardized MSE of this interpolation method is given by81 + 1A f x A f t m T ; x − t )dxdt 2 ∫t ∈U ∫x ∈Z { 0 ( )}{ 0 ( )} (
0MSE[ S]
2πAB2lim= lim µ2
Z 0 →∞
Z
→∞
0
B

− 1
{A f (x)}m(T ; x − t )dxdt
 πAB2 ∫t ∈U ∫x ∈Z 0 0
1

1 +

 B 2πA2 B2 ∫t ∈U ∫x ∈U m(T ; x − t )dxdt
+
µ
− 1
{A f (x)}m(T ; x − t )dxdt  πA2 B2 ∫t ∈U ∫x ∈U 0
= Q2 µ 2 + Q1µ(15)(for details, see Sadahiro 1999b). Similar to the point-in-polygon method, the accuracy of
the areal weighting interpolation is measured by Q1 and Q2.
4.1 Size of the target zone
We first analyze the relationship between the size of the target zone and the accuracyof estimates. Intuitively, it is expected that the relative accuracy increases as T becomeslarge (Langford et al. 1991). We test this hypothesis on a square lattice using theinterpolation models.
For the target zone T, we adopt circles whose radii range from 0.1 to 5.0. The areaof square cells is set to 1. Three forms of f(x) are investigated: a uniform distribution(Figure 7a), a quadrangular pyramid distribution (Figure 7b), and a quadrangular prism
distribution (Figure 7c). The latter two are adopted as approximations of points located on
the representative point, the point distribution which the point-in-polygon method
implicitly assumes.
FIG. 7. Forms of f(x) in a cell. a) a uniform distribution, b) a quadrangular pyramid
distribution, c) a quadrangular prism distribution.
The results are shown in Figures 8 and 9. As we have expected, the relative
accuracy of estimates improves as T becomes larger. This tendency is found in bothinterpolation methods for any form of f(x). We should note, however, that Q2 changes
periodically with the size of T, whereas Q1 decreases monotonously with as increase of
T. This is mainly because of the periodicity lying in both f(x) and the lattice. This implies
that the MSE may change drastically with a slight enlargement or reduction of T if points
are regularly distributed like a population distribution in new towns.
FIG. 8. Estimation accuracy of the point-in-polygon method and the size of the target zone

9T.
FIG. 9. Estimation accuracy of the areal weighting interpolation method and the size of thetarget zone T. The value of Q2 is constant at zero for the uniform distribution of f(x).
4.2 Distribution of points - the degree of concentration
We next investigate the relationship between the distribution of points and theaccuracy of estimates. The point-in-polygon method implicitly assumes that points arelocated exactly on representative points as mentioned earlier. The areal weighting
interpolation method, on the other hand, assumes that points are uniformly distributed inzones. We expect that an interpolation method yields accurate estimates if a point
distribution is close to its underlying assumption. To confirm this we analyze how theconcentration of points affects the accuracy of the point-in-polygon method.
For the target zone T we try three sets of circles whose radii range from 0.62-1.12,
1.62-2.12, and 3.62-4.12 since the MSE changes periodically with the size of T as seenin the previous subsection. We calculate the averages of Q1 and Q2 for each set of T. The
cell size is again set to 1. The form of f(x) changes from the uniform distribution to highly
concentrated distribution as shown in Figure 10. The width of a quadrangular prism
representing f(x) is denoted by w.
FIG. 10. Form of f(x) in a cell.
The results are shown in Figures 11a, 11b, and 11c. We notice that the figures arevery similar though the vertical scale is different. This indicates that the relationshipbetween the concentration of points and estimation accuracy is independent of the size ofthe target zone. Concerning the point-in-polygon method, Q1 and Q2 decreases as points
concentrate around the representative point. This supports the hypothesis mentioned
earlier. In contrast to this, the areal weighting interpolation method is not fully consistent
with the hypothesis: as w decreases, Q2 increases but Q1 decreases. This implies that if
the point density µ is low enough (say, µ=1) estimation accuracy improves as the point
distribution becomes dissimilar to the underlying distribution of the areal weighting
interpolation method. Though such a low density of points is not frequently observed inpractice, it should be noted that estimation error does not always decrease when a point
distribution becomes close to the underlying assumption of the interpolation method.
We now consider the degree of concentration of points where the two interpolation
methods give the same accuracy. Interestingly, the two methods are equivalent inestimation accuracy when w=0.35-0.38 for both Q1 and Q2 for any size of T. From this10

we can say that the point-in-polygon method is more suitable than the areal weighting
interpolation method when points are concentrated in less than 12-15% area around therepresentative point. If points are supposed to be more dispersed, the areal weighting
interpolation method would be appropriate.
FIG. 11. Estimation accuracy and the concentration of points. Radius of T is a) 0.621.12, b) 1.62-2.12, c) 3.62-4.12.
4.3 Distribution of points - relative location to the representative point
As assumed in this section, the representative point is often located on the centroidof a source zone, especially in digital spatial data. The location of points, however, does
not usually agree with representative points: points are sometimes dispersed in a zone,
clustered apart from a representative point, or even located along the edge of a zone. Sucha locational gap between points and the representative point is a source of estimation error
in the point-in-polygon interpolation. Therefore, we finally analyze how the location ofpoints in cells affects estimation accuracy.
To this end, we employ three quadrangular prism distributions as f(x) whose
widths are 0.1, 0.3, and 0.5. Points move gradually from the centroid of cells to thecorner as shown in Figure 12. The distance between the center of prism and therepresentative point located on the centroid of a cell is denoted by d. For the target zone T
we try the same set of circles as those used in the previous subsection, and calculate theaverages of Q1 and Q2. The area of cells is set to 1.
FIG. 12. Location of points in a cell.
As seen in Figure 13, the obtained results do not much differ according to the sizeof the target zone. In general, estimation error increases consistently with d. It isinteresting that this tendency is observed not only with the point-in-polygon method butwith the areal weighting interpolation method (Q1). The effect of locational gap is more
serious in the point-in-polygon method, especially when points are highly concentrated.
The point-in-polygon appears not so robust against the locational gap between points andthe representative point.
When points are located distantly from the representative point, it often happens thatthe areal weighting interpolation method gives better estimates even if points are tightly
clustered. The critical value of d where the two methods are equivalent in estimation
accuracy depends on w, the degree of concentration of points. Figure 13 depicts thatconcentration of points increases the critical value of d.

11

FIG. 13. Estimation accuracy and the location of points. Radius of T is a) 0.62-1.12, b)
1.62-2.12, c) 3.62-4.12.
5. OPTIMAL LOCATION OF REPRESENTATIVE POINTSIn the previous section we analyzed the effect of the locational gap between points
and the representative point. It was found that the point-in-polygon method is sensitive tothe locational gap, and that the estimation error increases with the distance. This implies
that if providers of spatial data such as national mapping agencies put the representative
point on an appropriate position it would improve estimation accuracy of the point-inpolygon interpolation performed by users of the data. So where should the representative
point be located in a zone? This section briefly discusses this subject.
We consider the optimal location of the representative point individually for eachzone, and thus we focus on the representative point of zone Zk. From the viewpoint ofareal interpolation, it is desirable that the representative point is located so as to minimize
estimation error originating in Zk. This problem is formulated as follows.
Let Mk be the number of points in the intersection of Zk and T. It is mathematically
represented as(16)
Mk = ∑ C y j U k y j .

( ) ( )j

The estimate of Mk given by the point-in-polygon interpolation isM̂k = C(pk )∑ Uk y j .

( )

j

(17)

The mean square error of Mk is calculated in a similar way as that of M (for details, see
Appendix 2):
MSEk [ S] = N ( N − 1)∫t ∈Z k∫x ∈Z kf ( x ) f (t ) Pr[ x ∪ t ∈ T ]dxdt Bf ( x )dx
+2 N   ∫ A0  x ∈Z k
−2 N ( N − 1)∫x ∈Z k−2 N ∫x ∈Z kf ( x )dx ∫x ∈Z kf ( x ) Pr[pk ∪ x ∈ T ]dx

f ( x ) Pr[pk ∪ x ∈ T ]dx

{

.

(18)

}

2
 B
f ( x )dx
+ N ( N − 1)  ∫ A0  x ∈Z i

The optimal location of pk is obtained by solving
min MSEk [ S].
pk

Substitution of equation (18) yields

12

(19)

min MSEk [ S] ⇔ max ∫pk

pk

x ∈Z kf ( x ) Pr[pk ∪ x ∈ T ]dx .

(20)Assuming that the region Z0 is even larger than T, we obtainmax ∫x ∈Z kpk

f ( x )m(T ; p k − x )dx .

(21)

The above problem is a kind of continuous location problems studied in operations
research (Plastria 1995). If f(x) and m(T; l) are given in analytical and simple forms, theproblem may be solvable by a computational procedure (Love 1972; Drezner andWesolowsky 1980; Aly and Marucheck 1982; Love et al., 1988). Otherwise we have to
try numerous locations in Zk to find the optimal location of pk.
Equation (21) implies that the optimal location of pk depends not only on f(x) andZk but also on the target zone T. In practice, however, the target zone is unknown when
the location of pk is determined. We hence consider the circle which is the most typicalshape of the target zone.
Let T be a circle of radius r whose diameter 2r is larger than the maximum length of
Zk. The measure m(T; l) is then given bylm(T ; l ) = 4πr 2 arccos  − πl 4r 2 − l 2 .
(22)
 2r 
Applying a Taylor series expansion around l=0 to equation (22), we obtaind1 d2m(T ; l ) ≈ m(T ; 0) + m(T ; 0)l +m(T ; 0)l 2
2
.
dl
2 dl

(23)

= 2π r − 4πrl
The above equation is a good approximation of m(T; l) especially for small l (Figure 14).
Substitution of equation (23) into equation (21) yieldsmin ∫f ( x ) p k − x dx .
(24)
2 2

pk

x ∈Z kThis indicates that the location of pk is the spatial median of f(x) in Zk (Brown 1983;
Small 1990; Stoyan and Stoyan 1994). As we will illustrate later, the spatial median does
not always agree with the centroid of f(x) on which the representative point is often
located.
FIG. 14. The measure m(T; l) of a circle of radius 3.
The above discussion holds for the shapes that allow the linear approximation
shown in equation (23). Convex shapes similar to the circle, say, the square, regulartriangle, regular hexagon, and rectangles meet this requirement (see Figures 15 and 16).
Considering that the target zone usually has a simple convex shape, we may say that the
spatial median is the optimal location of representative points in relation to estimation
accuracy of the point-in-polygon interpolation.

13

FIG. 15. The measure m(T; l) of a regular triangle and regular hexagon.
FIG. 16. The measure m(T; l) of rectangles having a variety of horizontal to verticalratios.
We finally show some examples of the optimal location of representative points.
We compute both the spatial median and the centroid of f(x) in square zones varying the
distribution of points f(x). Every zone is divided into 10 * 10 cells each of which has aconstant value of f(x) as shown in Figure 17.
As mentioned earlier, the spatial median does not always agree with the centroid.
The spatial median is more strongly drawn by points, and it tends to be nearer the
maximum of f(x). Since the locational gap between points and the representative point is
crucial in the point-in-polygon interpolation, estimation accuracy appears to be fairly
improved if data providers locate the representative point on the spatial median rather thanthe centroid of f(x).
FIG. 17. The optimal location of the representative point (spatial median) in a squarezone. The white squares and black circles represent the spatial median and the centroid off(x), respectively. The value of f(x) is zero in the white cells. The lightest gray cells
indicate f(x)=1.0, and the values of f(x) in the other gray cells are shown individually.
6. CONCLUSION
In this paper we have analyzed the accuracy of count data estimated by the point-inpolygon method. We first proposed a point-in-polygon interpolation model which isbased on a stochastic distribution of points and the target zone. This allows us to replace
diverse situations by their representative model. We then applied the model to the case ofa square lattice and numerically investigated the accuracy of estimates. The major results
are summarized as follows:
1) The relative accuracy of estimates generally increases with the size of the target
zone. However, the periodicity in the spatial configuration of source zones andf(x) disturbs monotoneity of the relationship, thus an enlargement of the target
zone may reduce estimation accuracy.
2) The point-in-polygon method gives good estimates when points are concentrated
around the representative point. The areal weighting interpolation method, on theother hand, fits uniform distributions of points. These methods are equivalent interms of estimation accuracy when points are concentrated in less than 12-15%

14

area around the representative point.
3) The point-in-polygon method is not so robust against the locational gap between
points and the representative point. Hence it often yields worse estimates than
the areal weighting interpolation method even if points are strongly concentrated.
The third result implies that the location of representative points is crucial for the
point-in-polygon interpolation. Thus we discussed the optimal location of representative
points and found that in most cases it is given by the spatial median of points. This warns
us as data producers that the centroid of points often used as the representative point is
not the optimal location in relation to the point-in-polygon interpolation.
In the empirical study we considered only the square lattice due to limitations ofspace. We should emphasize, however, that the proposed method is applicable to any
zonal system as mentioned earlier. Zonal systems which are less artificial than regularlattices, say, census tracts and municipal districts, should be analyzed in future
researches.
The point-in-polygon interpolation model we proposed considers data transfer frommultiple source zones to a single target zone. Data transfer between zonal systems still
remains to be analyzed where estimation occurs in multiple target zones. This is an
important subject in geography which is related to the modifiable areal unit problem(Openshaw 1984; Flowerdew and Green 1991; Fisher and Langford 1995; Okabe andSadahiro 1997). We will examine an extension of our study to this case in order toanalyze the accuracy of estimation between zonal systems.

15

LITERATURE CITED
Aly, A. A. and A. S. Marucheck (1982). "Generalized Weber Problem with Rectangular
Regions." Journal of the Operational Research Society 33, 983-989.
Brown, B. M. (1983). "Statistical uses of the Spatial Median." Journal of the Royal
Statistical Society, Series B 45, 25-30.
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.
Drezner, Z. and G. O. Wesolowsky (1980). "Optimal Location of a Demand Facility
Relative to Area Demand." Naval Research Logistics Quarterly 27, 199-206.
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.
Flowerdew, R. (1988). "Statistical Methods for Areal Interpolation: Predicting Count
Data from a Binary Variable." Research Report 15, Northern Regional ResearchLaboratory.
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.
Love, R. F. (1972). "A Computational Procedure for Optimally Locating a Facility with
Respect to Several Rectangular Regions." Journal of Regional Science 12, 233-242.
Love, R. F., J. G. Morris and G. O. Wesolowsky (1988). Facilities Location: Models& Methods. Amsterdam: North-Holland.
Markoff, J. and G. Shapiro (1973). "The Linkage of Data Describing Overlapping

16

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.
Openshaw, S. (1984). The Modifiable Areal Unit Problem. Norwich: Geobooks.
Plastria, F. (1995). "Continuous Location Problems." In Facility Location: A Survey ofApplications and Methods, edited by Z. Drezner, pp.225-262, New York: SpringerVerlag.
Preparata, F. P. and M. I. Shamos (1985). Computational Geometry - An Introduction-. New York: Springer-Verlag.
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.
Ripley, B. D. (1981). Spatial Statistics. New York: John Wiley.
Sadahiro, Y. (1999a). "Statistical Methods for Analyzing the Distribution of Spatial
Objects in Relation to a Surface." 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:
Addison-Wesley.
Sinfonica (1994). Study on Mesh Data. Report.
Small, C. G. (1990). "A Survey of Multidimensional Medians." International Statistical
Review 58, 263-277.
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.

17

APPENDIX 1
The mean square error of ε which is denoted by MSE[S] is calculated as follows.
MSE[ S] = E ε 2.
(A 1)
ˆ +E M
ˆ2= E M 2 − 2 E MM

[ ]
[ ]

[ ] [ ]

The first term of equation (A 2) becomes
2

 
2
E M = E  ∑ C y j  
 j
 

( )

[ ]



= E ∑ ∑ C y j C y j '  j j'
( )( ).

[ ( ) ( )](A 2)

[ ( )]= 2 ∑ E C y j C y j' + ∑ E C y j
j ≠ j'

Substituting
E C y j C y j'j'

[ ( ) ( )] = Pr[(y ∈T ) ∩ (y ∈T )]= ∫ ∫ Pr[(y ∈ dx ) ∩ (dx ∈ T ) ∩ (y ∈ dt ) ∩ (dt ∈ T )]dxdtjt ∈Z 0=∫t ∈Z 0j'j

x ∈Z 0∫x ∈Z 0j'f ( x ) f (t ) Pr[ x ∪ t ∈ T ]dxdt
(A 3)

and

[ ( )] [E C y j = Pr y j ∈ T=∫

x ∈Z 0=]f ( x ) Pr[ x ∈ T ]dx ,

(A 4)

B
A0

we have[ ]E M 2 = N ( N − 1)∫t ∈Z 0∫x ∈Z 0 B
f ( x ) f (t ) Pr[ x ∪ t ∈ T ]dxdt + N   .
 A0 

(A 5)

The second term of equation (A 1) is rewritten asˆ = E  C y  C(p ) U y 
E MM∑ j ∑i i ∑j i j
 j



[ ]

( )

( )



= E ∑ ∑ ∑ C(pi )Ui y j C y j' 
 i j j'( ) ( )

[ ( )] [

( )]

.

[ ( )

( )]

= 2 ∑ ∑ E Ui y j E C(pi )C y j' + ∑ ∑ E Ui y j C(pi )C y jij ≠ j'

Substitution of

18

i

j

(A 6)

[ ( )] [( )] = ∫f ( x )dx ∫Pr (pi ∈ T ) ∩ y j' ∈ dx ∩ (dx ∈ T ) dx=∫

f ( x )dx ∫f ( x ) Pr[pi ∪ x ∈ T ]dx

E Ui y j E C(pi )C y j'x ∈Z ix ∈Z ix ∈Z 0x ∈Z 0

[

(

]

)

(A 7)
and

[ ( )

( )]

[(

x ∈Z iPr y j ∈ dx ∩ (pi ∈ T ) ∩ (dx ∈ T ) dxjx ∈Z i=∫x ∈Z iyields[ ]ˆ = N ( N − 1)
E MM∑∫

x ∈Z ii+N∑ ∫x ∈Z ii])= ∫ Pr[y ∈ dx ] Pr[p ∪ dx ∈ T ]dx

E Ui y j C(pi )C y j = ∫

(A 8)

i

f ( x ) Pr[pi ∪ x ∈ T ]dx
f ( x ) Pr[pi ∪ x ∈ T ]dx

f ( x )dx ∫x ∈Z 0.

f ( x ) Pr[pi ∪ x ∈ T ]dx(A 9)

The third term of equation (A 1) is given by
2



2
ˆ = E   C( p ) U y  EM∑ i ∑j i j 
 i
 

[ ]

( )



= E ∑ ∑ C(pi )C(pi' )∑ Ui y j ∑ Ui' y j' jj'
 i i'


( )

[

]

( )

.

(A 10)[ ( ) ( )]= ∑ ∑ E C(pi )C(pi' ) ∑ ∑ E Ui y j Ui' y j'ii'This equation becomes
ˆ2 =2
EM

[ ]

j

j'

∑ ∑ E[C(p )C(p )]∑ E[U (y )]E[U (y )]
i

i

i'

i

j ≠ j'

i'

[

]

j

i'

j'

[ ( ) ( )]+ ∑ ∑ E C(pi )C(pi' ) ∑ E Ui y j Ui' y j
i

i'

j

= N ( N − 1)∑ ∑ Pr[pi ∪ pi' ∈ T ]∫
i

[

i'

]

x ∈Z if ( x )dx ∫f ( x )dxf ( x )dx ∫f ( x )dxx ∈Z i '[ ( )]+ ∑ E C(pi ) ∑ E Ui y jij= N ( N − 1)∑ ∑ Pr[pi ∪ pi' ∈ T ]∫i

x ∈Z ii'

 B
+ N 
 A0 
Using equations (A 5), (A 9), and (A 11), we obtain

19

x ∈Z i '.

(A 11)

MSE[ S] = N ( N − 1)∫t ∈Z 0∫x ∈Z 0f ( x ) f (t ) Pr[ x ∪ t ∈ T ]dxdt

 B
+2 N   A0 
−2 N ( N − 1)∑ ∫x ∈Z ii−2 N ∑ ∫i

x ∈Z if ( x )dx ∫x ∈Z 0f ( x ) Pr[pi ∪ x ∈ T ]dx.

f ( x ) Pr[pi ∪ x ∈ T ]dx+ N ( N − 1)∑ ∑ Pr[pi ∪ pi' ∈ T ]∫i

x ∈Z ii'f ( x )dx ∫f ( x )dxx ∈Z i 'APPENDIX 2
The mean square error of Mk is calculated as follows.
ˆ +E Mˆ2
MSEk [ S] = E Mk2 − 2 E Mk M
k
k

[

[ ]

] [ ]

The first term of equation (A 13) becomes
2

 
2
E M k = E  ∑ C y j U k y j  
 j
 

(A 12)(A 13)( ) ( )[ ][ ( ) ( ) ( ) ( )][ ( ) ( )].

(A 14)= 2 ∑ E C y j Uk y j C y j ' Uk y j ' + ∑ E C y j Uk y jj ≠ j'Usingj[ ( ) ( ) ( ) ( )]= ∫ ∫ Pr[(y ∈ dx ) ∩ (dx ∈ T ) ∩ (y ∈ dt ) ∩ (dt ∈ T )]dxdt
E C y j Uk y j C y j ' Uk y j 't ∈Z k=∫t ∈Z kandx ∈Z k∫x ∈Z kjj'f ( x ) f (t ) Pr[ x ∪ t ∈ T ]dxdt[ ( ) ( )]E C y j Uk y j = ∫x ∈Z k=∫x ∈Z k=we have[ ](A 15)E Mk2 = N ( N − 1)∫

[(

)

]

Pr y j ∈ dx ∩ (dx ∈ T ) dx
f ( x ) Pr[ x ∈ T ]dx,(A 16)Bf ( x )dx
A0 ∫x ∈Z kt ∈Z k∫x ∈Z kf ( x ) f (t ) Pr[ x ∪ t ∈ T ]dxdt

 B+ N  ∫f ( x )dx
 A0  x ∈Z k
The second term of equation (A 13) is rewritten as

20.

(A 17)

ˆ = E C(p )E Mk Mkk ∑ ∑ C y j Uk y j Uk y j ' 
j
j'



[

]

( ) ( ) ( )= 2 E  C ( p k ) ∑ C y j U k y j U k y j '  + E  C ( p k )∑ C y j U k y j  .
j ≠ j'j( ) ( ) ( )[( ) ( )( ) ( )] [ ( )][( ) ( )]

= 2 ∑ E C(p k )C y j Uk y j E Uk y j' + ∑ E C(pk )C y j Uk y jj ≠ j'j(A 18)
Substituting equation (A 8), we obtainˆ = N ( N − 1)
E Mk Mk∫ f (x)dx ∫[]x ∈Z k+N∫

x ∈Z kf ( x ) Pr[p k ∪ x ∈ T ]dx

x ∈Z k.

f ( x ) Pr[p k ∪ x ∈ T ]dx(A 19)The third term of equation (A 13) becomes
2

 
2
ˆ
E Mk = E C(pk )∑ Uk y j  

j
 

[ ]

( )





= 2 E  C ( p k ) ∑ U k y j U k y j '  + E  C ( p k )∑ U k y j j ≠ j'
j





[

]

( ) ( )( )[ ( )] [ ( )] [][ ( )]= 2 E C( p k ) ∑ E Uk y j E Uk y j ' + E C( p k ) ∑ E Uk y jj ≠ j'{= N ( N − 1) Pr[pi ∈ T ] ∫{}f ( x )dxx ∈Z i}2j+ N Pr[pi ∈ T ]∫x ∈Z if ( x )dx

. (A 20)

2
 B
 B
f ( x )dx + N   ∫= N ( N − 1)  ∫f ( x )dx
 A0  x ∈Z i
 A0  x ∈Z i

Using equations (A 17), (A 19), and (A 20) we obtain
MSEk [ S] = N ( N − 1)∫t ∈Z k∫x ∈Z kf ( x ) f (t ) Pr[ x ∪ t ∈ T ]dxdt Bf ( x )dx
+2 N   ∫ A0  x ∈Z k
−2 N ( N − 1)∫x ∈Z k−2 N ∫x ∈Z kf ( x )dx ∫x ∈Z kf ( x ) Pr[pk ∪ x ∈ T ]dx

f ( x ) Pr[pk ∪ x ∈ T ]dx

{

}

2
 B
f ( x )dx
+ N ( N − 1)  ∫ A0  x ∈Z i

21.

(A 21)

Z0
Z1
Z3

a.

Z2

Z5

Z4

Z0
b.

(13)

p1

(10)

c.
(4)

Z0

p2

p3

(6) p
4

Figure 1

(3)

p5

(13)

(10)

a.
(4)

(3)
(6)
4 + 10 + 6 = 20

(13)

(10)

b.
(4)

(3)
(6)
4×

1
1
1
+ 10 ×
+6×
=5
2
10
3

Figure 2

Z0

Figure 3

T
Z0

Z0
T

Figure 4

t

Figure 5

Figure 6

4.0
3.0

1.0

1.0

a.

0.5

1.0

b.

Figure 7

c.

102
Q1

101

Q2

f(x)
uniform
quadrangular pyramid
quadrangular prism

100
10-1
Q1 and Q2

10-2
10-3
10-4
10-5
10-6
0.0

1.0

2.0

3.0
radius of T

Figure 8

4.0

5.0

102
Q1

101

Q2

f(x)

(zero)

uniform
quadrangular pyramid
quadrangular prism

100
10-1
Q1 and Q2

10-2
10-3
10-4
10-5
10-6
0.0

1.0

2.0

3.0
radius of T

Figure 9

4.0

5.0

1.0/w2w

Figure 10

0.3
Q1

Q2

interpolation method
point-in-polygon
areal weighting

0.2

Q1 and Q2
0.1

0.0
0.0

0.2
0.4
0.6
0.8
width of the quadrangular prism distribution (w)Figure 11a

1.0

0.03
Q1

Q2

interpolation method
point-in-polygon
areal weighting

0.02

Q1 and Q2
0.01

0.00
0.0

0.2
0.4
0.6
0.8
width of the quadrangular prism distribution (w)Figure 11b

1.0

0.003Q1

Q2interpolation method
point-in-polygon
areal weighting

0.002

Q1 and Q2
0.001

0.000
0.0

0.2
0.4
0.6
0.8
width of the quadrangular prism distribution (w)Figure 11c

1.0

w

Figure 12

0.5
Q1

Q2

0.4

interpolation method
point-in-polygon
areal weighting

0.3
Q1 and Q2
0.2

0.1

0.0
0.00

0.10

0.20

0.30

0.40 0.45 0.00

0.10

0.20

0.30 0.35 0.00

0.10

distance between the center of f(x) and the representative point (d)

w=0.1

w=0.3

Figure 13a

w=0.5

0.20 0.25

0.05
Q1

Q2

0.04

interpolation method
point-in-polygon
areal weighting

0.03
Q1 and Q2
0.02

0.01

0.00
0.00

0.10

0.20

0.30

0.40 0.45 0.00

0.10

0.20

0.30 0.35 0.00

0.10

distance between the center of f(x) and the representative point (d)

w=0.1

w=0.3

Figure 13b

w=0.5

0.20 0.25

0.005
Q1

Q2

0.004

interpolation method
point-in-polygon
areal weighting

0.003
Q1 and Q2
0.002

0.001

0.000
0.00

0.10

0.20

0.30

0.40 0.45 0.00

0.10

0.20

0.30 0.35 0.00

0.10

distance between the center of f(x) and the representative point (d)

w=0.1

w=0.3

Figure 13c

w=0.5

0.20 0.25

7.0

6.0

5.0

4.0
m(T; l)
3.0

2.0

1.0

0.0
0.0

2.0

4.0
l

Figure 14

6.0

8.0

7.0

6.0

5.0

4.0
m(T; l)regular hexagon
3.0

2.0

regular triangle

1.0

0.0
0.0

2.0

4.0
l

Figure 15

6.0

8.0

7.0

6.0

5.0

4.0
m(T; l)

1:1
3.0

1:3

2.0
1:5
1.0

0.0
0.0

2.0

4.0
l

Figure 16

6.0

8.0

1.1

1.2

1.5

2.0

3.0

1.1

1.2

1.5

2.0

3.0

Figure 17