Introduction
Monitoring the status of the global environment
and assessing the condition of ecological resources across the earth
have emerged as major scientific activities for the new century.
Environmental phenomena and problems operating at different spatial
and temporal scales, biogeochemical cycles and climate change being
but two examples, are now being studied as global systems. A key
issue is how to integrate both spatially and temporally disparate
global data now being collected and made publicly available to the
scientific community (Brooks 1981; Kahn and Braverman 1999).
One approach to data integration is to partition
the earth into sampling or analysis units that form a hierarchy.
Partitioning based on a hierarchy of political divisions, such as
nations, states and counties, is a familiar approach linked to governmental
data collection and dissemination efforts. However, scientists often
favor developing a hierarchical, geometrically regular global partitioning
system that is unbiased with respect to spatial patterns created
by natural and human processes (White et al. 1998). A "flat
earth" example would be to divide the flat surface into identical
grid squares, each of which is then partitioned into quarters. This
is called a four-fold partition (Figure 1 left), whereas
partitioning an initial triangle into nine identical triangles produces
a nine-fold partition (Figure 1 right). The grid cell hierarchy
is formed by successive partitionings that we call recursion levels.
The geometrically ideal global partitioning system
would consist of grid cells equal in surface area and identical
in shape, akin to the square, equilateral triangular, or regular
hexagonal grid cells that tessellate a flat surface. The centerpoints
of these cells form geometrically uniform rectangular and triangular
point sampling grids advantageous in statistical sampling designs
for local areas. On the entire globe, however, area and shape regularity
can be achieved only by projecting the faces of a Platonic polyhedron
(tetrahedron, hexahedron, octahedron, dodecahedron, or icosahedron)
onto the spherical or ellipsoidal approximations to the earth's
true shape. Since further partitioning of a spherical or ellipsoidal
polyhedron face introduces unavoidable variation in cell area, shape,
or both, grid cell optimization with respect to area, shape, and
other criteria is a fundamental issue (Kimerling et al. 1999).
|
Figure 1. Tobler-Chen
cylindrical equal area and Icosahedral Snyder equal area map
projection surfaces, each covering 1/20th of the earth's surface.
The Tobler-Chen surface has been four-fold partitioned to
recursion level 3 and the Snyder surface has been nine-fold
partitioned to recursion level 3.
|
Numerous approaches to partitioning the globe
into areal grid cells have been proposed, and both Dutton (1999)
and Kimerling et al. (1999) have grouped partitioning schemes into
a small number of general categories. At the broadest level, partitioning
methods can be divided into a) those that directly subdivide the
sphere or ellipsoid into grid cells and b) those that project the
globe onto one or more map projection surfaces which are partitioned
in a geometrically regular manner and then back-projected to the
globe.
Equal area grid cells are a primary requirement
to many scientists and spatial statisticians. Of the map projection
based partitioning methods, the equal area polyhedral projections
derived by Snyder (1992) satisfy this need (Figure 1 right), as
does the cylindrical equal area projection approach (Figure 1 left)
of Tobler and Chen (1986). Both introduce unavoidable variation
in grid cell shape and compactness, the latter method suffering
far greater in this respect. Directly partitioning the sphere into
equal area quadrilateral cells has been achieved by mathematicians
(Rakhnanov et al.1994) for a small number of cell densities, at
the cost of inherently wide variations in cell shape characteristic
of "constant area" global partitions (Kahn and Braverman
1999). Similar geometrical irregularities in point spacing have
been noted by mathematicians studying the topologically dual problem
of packing points on a sphere optimally (Saff and Kuijlaars 1997:
Coexeter 1962). Directly partitioning the faces of Platonic polyhedra
cannot achieve the equal area criterion when each face is divided
into equal area sub-triangles and the three edges of each triangle
are geodesic (great circle) lines defined by the laws of spherical
trigonometry (Baumgardner and Frederickson. 1985).
To better satisfy the grid evaluation criteria,
a new global partitioning method, the small circle subdivision
method, is proposed. The small circle subdivision method, based
upon the direct spherical subdivision approach, partitions the sphere
or ellipsoid with small circles into sub-triangles equal in area
that vary only slightly in shape. We first examine in detail the
mathematical basis of the small circle subdivision method applied
to a spherical icosahedron face, and then present results from implementing
the method over a wide range of cell sizes. Results are presented
in the form of a comparison with another proposed equal area global
partitioning method based on the Icosahedral Snyder Equal Area (ISEA)
map projection.
[top]
Small Circle Subdivision Method
The small circle subdivision method directly partitions
the surface of the earth with small circles into equal area sub-regions.
The small circle subdivision method can eliminate the inflexibility
of great circle subdivision, because only one great circle is available
between two points on a sphere (that is, the curvature of the great
circle is fixed). Small circles can be any circles on a sphere with
an arbitrary radius and axis. The curvature of a small circle through
two points on a sphere can be flexibly changed. Thus, it is always
possible to partition a triangle into equal area sub-triangles if
appropriate small circles are computed.
Before describing the method of small-circle subdivision,
terms that will be used frequently in the section must be defined.
A circle on a sphere is an intersection of a sphere and a plane.
A small circle is defined as any circle on a sphere. If the
center of a small circle coincides with the center of the earth,
we call it a great circle. Obviously, great circles are a
subset of small circles. The axis of a small circle is the
straight line that passes through the earth's center and is perpendicular
to the plane of the small circle. The pole of a small circle
is the point at which its axis intersects the surface of the earth
(there are two intersection points, but we only consider the point
close to the small circle plane as its pole). There exists a great
circle parallel to the plane of a small circle, which we call the
relative equator of the small circle. We define the relative
latitude of a small circle as the angle between the plane of
its relative equator and the straight line which passes through
the sphere center and any point on the small circle. The relative
co-latitude of a small circle is the angle between the axis
of the small circle and the straight line through the sphere center
and any point on the small circle. Figure 2 shows a small
circle cbd, its axis AA¢,
its pole P, its relative equator CBD, its relative latitude ÐBOb
and its relative co-latitude ÐbOP.
Obviously, the axis, pole, equator and geodetic latitude in the
traditional spherical coordinate system are particular examples
of the above definitions. We also define
|
Figure 2. A small circle on a
sphere with axis AA¢,
pole P, relative equator CBD, relative latitude ÐBOb
and relative co-latitude ÐbOP.
|
a small circle triangle on a sphere as
a region that is on the sphere and confined by three small circles
(Figure 3). If all three circles are great circles, the triangle
is called a spherical triangle. Obviously, spherical triangles
are a subset of small circle triangles.
|
Figure 3. A small circle triangle on
a sphere.
|
Our goal is to partition a particular spherical
triangle, an icosahedron face, into equal area sub-triangles.
The edges of the sub-triangles could be either small or great circles.
We will consider two ways of dividing a triangle: four-fold and
nine-fold subdivision. In the four-fold subdivision method, the
three edges of a triangle are bisected and the bisection points
are connected with small circles, producing four sub-triangles within
the triangle. In the nine-fold subdivision method, the three edges
of a triangle are trisected , and the trisection points as well
as an inner center point are connected with one another by small
circles, producing nine sub-triangles within the triangle. We define
the initial spherical triangle (e.g. the icosahedron face) as a
recursion level 0 triangle. If a recursion level 0 triangle
is partitioned, we define the process as a recursion level 1 process
and the resulting sub-triangles as recursion level 1 triangles.
If the recursion level 1 triangles are partitioned, we define the
process as a recursion level 2 process and the resulting sub-triangles
as recursion level 2 triangles, and so on.
[top]
Four-Fold Subdivision Method
In the four-fold subdivision method, the three
edges of a triangle are bisected and the bisection points are connected
with small circles, producing four sub-triangles within the parent
triangle. If the bisection points were connected with great circles,
the four resulting sub-triangles would not have equal areas because
only one great circle passing through two points on a sphere is
available, that is, the curvature of the great circle is fixed.
For example, if we partition an icosahedron face, the great circle
segment between two bisection points always bends outward from the
center of the parent triangle, resulting in a larger area of the
middle sub-triangle (Figure 4). If the bisection points are
connected with small circles, we could obtain four sub-triangles
equal in area, because there exist infinite numbers of small circles
between two points on a sphere and we can arbitrarily change the
radius and axis (or curvature) of a small circle by moving its pole
up or down on the sphere. In Figure 4, if we move the pole
of a small circle toward the center sub-triangle, the radius of
the small circle will become smaller and the arc of the small circle
segment will shift toward the center of the parent triangle, resulting
in a smaller area for the center sub-triangle and a larger area
for the outer sub-triangle associated with the small circle segment.
The area of the outer sub-triangle is determined only by its associated
small circle segment, independent of the other two small circle
segments. If the pole of the small circle is moved to a position
on the sphere such that the area of the outer sub-triangle is just
one-fourth of the area of its parent triangle, the small circle
segment is determined. The other two small circle segments can be
determined in the same way. Because each of the four sub-triangles
has an area one-fourth of the parent triangle, they obviously are
equal in area.
|
Figure 4.
Four-fold small circle and great circle subdivision geometry
for an icosahedron face.
|
[top]
Mathematical Derivations
We now derive mathematical equations and determine
numerical methods to implement the subdivision process discussed
above, so that the following essential information about sub-triangles
at each recursion level will be obtained: (1) the latitudes and
longitudes of all vertices; (2) the relative co-latitudes of all
edges; (3) the geodetic latitudes and longitudes of poles for all
edges; (4) the arc lengths of all edges; and (5) the area of all
sub-triangles.
For simplicity, we represent the earth as a unit
sphere, so that the arc length (radian measure) of a great circle
segment between two points on the sphere is the same as the angle
it subtends. For a sphere with radius r, linear measurements are
scaled by a factor of r and area measurements by r2.
|
Figure 5. Four-fold small circle subdivision
for an arbitrary triangle.
|
Consider an arbitrary initial triangle ABC on
a unit sphere in Figure 5, with its three edges made of either
small or great circles. We assume that the following information
about the triangle is known a priori:
(1) the latitude and longitude of each
vertex;
(2) the latitude and longitude of the pole for each side;
(3) the arc length of each side (small circle segment);
(4) the relative co-latitude of each side.
(5) the area of the triangle.
Our problem is to find three small circles connecting
each pair of the midpoints D, E and F, so that four resulting sub-triangles
will have equal areas. We first find the coordinates of midpoints
D, E and F. Consider, in particular, the edge AB, with its pole
PAB at 
and its relative co-latitude aAB
(Figure 6). Denote the spherical coordinates for the vertices
A and B as 
and  ,
respectively. Construct a Cartesian coordinate system with the origin
located at the center of the earth, the x-y plane on the equator
and the z-axis pointing to the north pole. Standard formulas for
transforming from spherical coordinates to Cartesian coordinates
of a point on the unit sphere are:
(1)
The Cartesian coordinates,  ,
 and
 for
points A, B, and PAB, respectively, can be computed using
equations (1). Thus, the Cartesian coordinates of the midpoint
D¢
of the chord between the point A and the point B are:
(2)
|
Figure 6.
Geometry of a small circle midpoint calculation.
|
We denote the vector from the earth center O to
the point O¢,
the center of the small circle, as  ,
and the vector from the earth center O to the pole PAB
of the small circle as  .
So we get:
(3)
or
(4)
The vector 
which points from the earth center O to the midpoint D of the small
circle segment AB is given by:
(5)
where 
is a unit vector pointing from the center O¢
of the small circle to the
point D¢
or D, and sinaAB
is the radius of the small circle. Transforming equation (5)
from vector expression to scalar expressions, we get the Cartesian
coordinates of the midpoint D:
(6)
The spherical coordinates of the midpoint D are
calculated using the following equations:
(7)
The other two midpoints, E and F, can be found
in the same manner.
After three midpoints have been found, we next
determine three unknown small circles DE, EF and DF in Figure
5. First consider the small circle DF and the sub-triangle ADF.
We will find the small circle DF so that the area of sub-triangle
ADF is one-fourth of the area of the parent triangle ABC. We have
shown that the curvature of a small circle depends both on the position
of its pole and on its relative co-latitude. But for a small circle
crossing two fixed points on a sphere (i.e., the geodetic coordinates
of the two points are known), its relative co-latitude will be the
only variable in determining the small circle because its pole can
be found from the relative co-latitude. To find the relative latitude
of the small circle, we connect two midpoints, D and F, with a great
circle (Figure 7). Great circle segment DGF and small circle
segment DSF bound the shaded area denoted by SL. We call it a semi-lune
(in spherical trigonometry a lune is a region on a sphere bounded
by two great circle segments). Because the small circle segments
AD and AF and the great circle segment DGF are known, the area of
the triangle ADGF confined by the three segments can be computed.
Thus, the area of the semi-lune SL, denoted by SSL
, is calculated as the area difference between the triangle ADGF
and the triangle ADSF that is one-fourth the area of the parent
triangle. To find the area of the triangle ADGF, construct a spherical
triangle AGDGFGA (Figure 8) by connecting each pair of the
points A, D and F with a great circle. The area of triangle ADGF
(i.e. the triangle ASDGFSA in Figure 8) is the area difference
or sum of the spherical triangle AGDGFGA and two semi-lunes (SL1
and SL2) depending on the curvature of the small circles
ASD and ASF. For convenience, we assume that the area of a semi-lune
has a positive value if its small circle edge is outside its great
circle edge,
|
Figure 7. A semi-lune confined by a
great circle and a small circle.
|
and a negative value otherwise. In the Figure
8, for example, the semi-lunes SL1 and SL2
have negative areas and the semi-lune SL has a positive area.
So the area of the semi-lune SL can be written as:
(8)
where Sabc
is the area of the parent triangle ABC, SADF
is the area of the spherical triangle ADF (i.e. the spherical
triangle AGDGFGA in Figure 8), 
and 
are areas of two semi-lunes SL1 and SL2,
respectively. Our problem simplifies to finding the relative co-latitude
of the small circle DSF
|
Figure 8.
Geometry of area representation of the unknown semi-lune
SL by areas of two known semi-lunes and one spherical triangle.
Heavy lines represent small circles and light lines represent
great circles.
|
so that the area of the semi-lune SL satisfies
equation (8). The following steps are needed to calculate the
area of the semi-lune, SSL
:
(1)
Compute the area, SADF,
of the spherical triangle ADF. The arc lengths of great circle
segments AGD, DGF and AGF denoted by f, a and d, respectively,
can be obtained by calculating the spherical distances between
pairs of vertices. The arc length f of great circle segment AGD,
for example, is:
or
(9)
where 
and 
are the latitudes of the vertices A and D, respectively, and 
and 
are the longitudes of the vertices A and D. Similarly, we can
obtain the arc lengths a and d of the great circle segments DGF
and AGF. Denote the vertex angles of the spherical triangle by
ÐA,
ÐD
and ÐF.
The law of cosines for sides of a spherical triangle gives:
(10)
So the vertex angle ÐA
is:
(11)
Similarly, the vertex angles, ÐD
and ÐF,
can be obtained. The area of the spherical triangle ADF is the
spherical excess (sum of the interior angles minus p),
or:
(12)
(2)
Compute the areas 
and 
of two known semi-lunes SL1 and SL2. First
consider semi-lune SL1 (Figure 9). The pole,
PAD, of the small circle ASD is located at 
and its relative co-latitude is aAD.
Connect points A and D with great circles to the pole PAD
and calculate the area of the spherical triangle PADAGD.
The area
|
Figure 9.
Geometry of area calculation for the known semi-lune SL1.
The
small circle AD has its pole at PAD
and relative co-latitude aAD.
|
calculation is similar to that above, i.e.,
first compute the vertex angles of the spherical triangle, ÐPAD,
ÐA
and ÐD,
and then calculate the area:
(13)
The area of the small circle triangle PADASD
(polar wedge) is given by:
(14)
Then the area of the semi-lune SL1
is:
or
(15)
Because the vertex angle ÐA
is equal to ÐD,
we get
(16)
Now we need to find the angles ÐA
and ÐPAD.
We know that the arc length d (or a) of the great circle segment
PADGA (or PADGD)
is the relative co-latitude, aAD,
of the small circle ASD. The arc length, pAD,
of the great circle segment AGD is the spherical distance from
point A to point D, which can be calculated using equation
(9). The law of cosines for sides of a spherical triangle
gives:
(17)
Because d=a=aAD,
equation (17) becomes:
(18)
So the vertex angle ÐPAD
is:
(19)
Now calculate the vertex angle ÐA.
The law of cosines for angles of a spherical triangle gives:
(20)
A manipulation of equation (20) gives
the vertex angle A:
(21)
Substituting equations (19) and (21)
into (16), the area of semi-lune SL1 is computed.
Similarly, we can compute the area of semi-lune SL2.
(3)
Find the area, SSL,
of the unknown semi-lune SL. After the area, SADF,
of the spherical triangle ADF and the areas, 
and  ,
of semi-lunes SL1 and SL2 are found, the
area of the semi-lune SL can be computed from equation (8).
Next, determine the small circle DSF in Figure
7. As mentioned above, if we find the relative co-latitude of
the small circle, the small circle will be determined. Denote the
relative co-latitude of the small circle as aDF
and its pole as PDF. Similar to
the area calculation above for known semi-lune SL1, the
area of the semi-lune SL can be expressed by aDF
as follows:
(22)
and
where pDF is
the spherical distance from the point D to the point F, which is
computed using equation (9). Angles ÐPDF
and ÐD
are vertex angles of the spherical triangle PDFDGF
(Figure 10). There are three equations and three unknowns (aDF,
ÐPDF
and ÐD),
so aDF
can be determined. It is very difficult to get an analytical solution
from these equations because the set of equations is nonlinear.
So we turn to numerical methods to get an approximate solution with
prescribed precision. We assume an initial value of aDF
and find the area of the semi-lune by solving equation (22).
Then we adjust the value of aDF
up or down by comparing the computed area with the desired amount
obtained above. If the computed area is larger or smaller than the
desired amount, the value of aDF
is increased or decreased, respectively. If the difference between
the computed area and the desired amount is less than a prescribed
precision, we stop computation and select the value of aDF
as the final relative co-latitude of
|
Figure 10. Geometry of relative
co-latitude calculation of the small circle DF
with its pole at PDF and
relative co-latitude aDF.
|
small circle DSF. To speed up the convergence,
we will apply the binary-tree method. We first find the range of
the relative co-latitude, and then determine the relative co-latitude
iteratively using bisection beginning with the middle point of the
range. Because we always use the pole which is close to the small
circle plane, the maximum value of the relative co-latitude is p/2.
If we are only interested in the small circle segment DSF less than
a semi-circle, the minimum value of the relative co-latitude is
arcsin(tDF/2) where tDF
is the chord length between D and F. This consideration is reasonable
because we require that sub-triangles have similar shapes. If one
of the small circle edges of a sub-triangle has an arc length more
than a semi-circle, the shape of the small circle triangle will
be unacceptably strange. This requirement sets a limitation on subdivision
precision. In some cases where a triangle has an unusual shape and
the prescribed precision is high, a small circle segment larger
than a semi-circle may be needed to subdivide the triangle. In such
cases, if we search the value of aDF
only in the range (arcsin(tDF/2),p/2),
the computation would loop forever and the small circle would never
be found. In practice, this situation has never been encountered
when partitioning an icosahedron face up to recursion level 9, but
was encountered once when partitioning an octahedron face at recursion
level 8.
After the relative co-latitude of the small circle DSF has been
found, we next determine the spherical coordinates of the pole of
the small circle,  ,
in Figure 11. Denote the Cartesian coordinate of the pole
PDF by  .
The pole PDF is on the sphere:
(23)
The Cartesian coordinates of the endpoints, D and F, denoted by
and  ,
respectively, can be obtained from equation (1). The angle
between the vectors 
and 
or 
is the relative co-latitude aDF
of the small circle DF (Figure 11), so the inner products
give:
(24)
Equations (23) and (24) are
solved together for  ,
and  .
There are two possible solutions here. If the small circle segment
DSF is outside the great circle segment
|
Figure 11. Geometry of pole determination
of the small circle DF.
|
DGF, we choose one solution close to the vertex
A. Otherwise, we choose the solution far from the vertex A. Then
we transform the Cartesian coordinates back to spherical coordinates:
(25)
The arc length of the small circle segment
DSF can also be obtained from its relative co-latitude and the vertex
angle ÐPDF.
The expression for the small circle arc length is:
(26)
We determined the small circle DF in
Figure 5. The other two small circles DE and EF can be obtained
in the same way as above. Each of the four sub-triangles has an
area equal to one-fourth the area of their parent triangle, so the
triangle ABC is divided into four equal-area sub-triangles.
[top]
Nine-fold Subdivision Method
In the nine-fold method, all three edges of a
triangle are trisected and the trisection points, as well as a central
inner point, are connected by small circles, producing nine sub-triangles
within the parent triangle. The basic problem is still to determine
the small circles so that the resulting nine sub-triangles have
equal areas. The method of small circle determination is similar
to the four-fold case, with a few exceptions. First, the central
inner point of a triangle needs to be found and to be connected
with six trisection points of three sides with small circles so
as to construct nine sub-triangles (Figure 12). If six trisection
points are directly connected with one another by small circles,
three inner small circles, each dividing the hexagon DEFGHI into
two equal-area quadrilaterals, may not intersect at one point. The
position of the inner point will affect the compactness of the subdivision
method. We have tried several ways to choose the inner point, and
will illustrate these in detail later. Secondly, each of six inner
sub-triangles has two unknown sides. This is different from the
two-frequency method where each sub-triangle has only one unknown
side. Thus, to determine six inner small
|
Figure 12. Small circles directly
connect six trisecting points with one another. Three inner
small circles (heavy lines), each of them dividing the hexagon
DEFGHI into two equal-area quadrilaterals, do not intersect
at one point and so an inner point is required.
|
circles in three-frequency subdivision method,
at least one initial small circle should be known a priori to simplify
computations. The others will be derived one by one from the initial
small circles.
Consider the small circle triangle ABC in
Figure 13. First, we need to find the trisection points for
the three sides of the triangle. Consider, in particular, the edge
AB, with its pole PAB at 
and its relative co-latitude aAB
(Figure 14). Denote the spherical coordinates for the vertices
A and B as 
and  ,
respectively. The Cartesian coordinates,  ,
and 
for points A, B, and PAB can be
|
Figure 13. Nine-fold small circle
subdivision for an arbitrary triangle ABC.
|
calculated from equation (1). Assume
that D and E are two trisection points of side AB. Because D and
E are on the unit sphere, their Cartesian coordinates 
and 
satisfy:
(27)
(28)
|
Figure 14. Geometry of trisecting
point calculation of a small circle AB with its pole at PAB
and relative co-latitude aAB.
|
We know that the angle between the vector 
and each of the vectors 
and 
is the relative co-latitude, aAB,
of the small circle AB, so the inner product gives:
(29)
(30)
Denote the angle between vectors 
and 
as d,
so the inner product gives:
(31)
The angle between vectors 
and 
is also d,
so we get:
(32)
We need to compute the angles d.
Applying the plane law of cosines to the plane triangle connecting
A, D and the earth center O (Figure 14), the chord length
between A and D is:
(33)
In the plane triangle connecting A, D and the
center O¢
of the small circle AD, the length O¢A
or O¢D
is the radius of the small circle AD, which is r=sinaAB
(Figure 14). The angle ÐAO¢D
is one-third of the angle ÐAO¢B
which is equal to the arc length of the small circle segment AB,
lAB, divided by r. So the plane
law of cosines gives:
or
(34)
Solving the equations (33) and (34)
for cosd:
(35)
Substituting into equations (31) and (32),
we get:
(36)
(37)
We must solve equations (27), (29) and
(36) for  ,
and equations (28), (30) and (37) for  .
Each has two solutions. We only select the solution close to the
point B as the Cartesian coordinates for point D, and the solution
close to the point A as the Cartesian coordinates for point E. The
trisection points of other two sides can be computed in the same
way.
After the trisection points for three sides of
triangle ABC are found, we next determine the three outer small
circles DI, FE and HG. The method to find the three small circles
has been introduced previously in the four-fold subdivision method.
The only difference is that each area of the three outer sub-triangles
ADI, BFE and CHG is one-ninth the area of triangle ABC.
Next, we need to determine the inner point of
triangle ABC. As mentioned previously, the determination of the
inner point will affect the compactness of the subdivision method.
So the position of the inner point M is dependent upon the shape
of the hexagon DEFGHI which is confined by six small circles. The
ideal way is to select the center point of the hexagon as the inner
point. But the calculation of the center point is complex because
it involves integration over the unit sphere bounded by six small
circles. The other way to locate an inner point is to find three
small circles EH, GD and IF using the same method as in the four-fold
case so that each of three triangles AEH, BGD and CIF has an area
four-ninths the area of triangle ABC. Actually, each of the three
small circles divides the hexagon DEFGHI into two equal-area quadrilaterals.
These three small circles intersect at three points and define a
small triangle (see Figure 12). We arbitrarily choose a point
inside the small triangle as an inner point (e.g., the coordinates
of the point could be the average coordinates of the three intersection
points). However, a difficulty will be encountered later when we
determine six inner small circles using the inner point. As we explain
later, at least one initial small circle should be known a priori
to find other small circles. But there is no proper way to obtain
such a small circle if we use the inner point. To simplify our calculation,
we use two other methods to find the inner point. First, find the
small circle EH that divides the hexagon DEFGHI into two equal-area
quadrilaterals, and then choose the bisection point of the small
circle segment EH as an inner point (Figure 15). Second,
find the small circles DG and IF in the same way as above, and then
select the intersection point of the two small circles as an inner
point (Figure 16). The Cartesian coordinates 
of the intersection point T, can be calculated by the following
steps:
(1)
Denote the relative co-latitude of small circle
DG as aDG,
the spherical coordinates of its pole PDG
as  ,
the relative co-latitude of small circle IF as aIF,
and the spherical
|
Figure 15. Nine-fold small circle
subdivision with the inner point being the midpoint of the
small circle EH which divides the small circle hexagon DEFGHI
into two equal-area quadrilaterals.
|
coordinates of its pole PIF
as  .
The Cartesian coordinates of the poles PDG
and PIF 
and  ,
can be calculated using equation (1).
(2)
The angle between the vector 
and the
vector 
is aDG
and the angle between
vector 
and 
is aIF.
So the inner product gives:
(38)
(39)
|
Figure 16. Inner point determination
for nine-fold small circle subdivision. Each of the small
circles DG and FI divides the hexagon DEFGHI into two equal-area
quadrilaterals.
|
The intersected point T is on the unit sphere,
so we get
(40)
Solving equations (38), (39)
and (40), we get two solutions for the Cartesian coordinates
of the intersection point T. We selected the solution close to vertex
A.
After the trisection points and the inner
point are found, we finally determine the six inner small circles
DT, ET, FT, GT, HT and IT. There are six variables for the six relative
co-latitudes of these small circles. If we apply the same method
as in the four-fold subdivision case to the six inner sub-triangles
so that each of them has an area one-ninth of the triangle ABC,
six equations can be obtained. But only five of the equations are
independent, so one more equation is needed to solve the six relative
co-latitudes. The sixth equation can be obtained by maximizing the
average compactness of the six sub-triangles. Although this method
is ideal and guarantees the equal area and maximum compactness,
it is practically impossible to implement because the six equations
are nonlinear. Instead, we assume at least one initial small circle
to be known a priori. The choice of the initial small circles depends
on the determination of the inner point. In Figure 15, if
we select the midpoint T of the small circle segment EH as the inner
point, we can assume the segments ET and HT to be two initial small
circles and their relative co-latitudes are the same as the relative
co-latitude of the small circle EH. Thus, the relative co-latitudes
of other four small circles can be easily calculated one-by-one
using the method introduced previously for four-fold subdivision.
In Figure 16, if we select the intersection point T of the
two small circles DG and IF as the inner point, one of the small
circle segments DT, GT, IT and FT can be selected to be an initial
small circle. The relative co-latitude of the initial small circle
is the same as that of the small circle DG or IF. In the figure,
we have selected the segment FT as an initial small circle. The
relative co-latitudes of other five small circles can be derived
one-by-one using the same method as in the four-fold subdivision
case.
After the relative co-latitudes of the nine small circles are computed,
the poles and the arc lengths of all the small circles can be calculated
the same way as in the four-fold case.
After the relative co-latitudes of the nine small
circles are computed, the poles and the arc lengths of all the small
circles can be calculated the same way as in the four-fold case.
|
Figure 17. Nine-fold small circle
subdivision with the initial small circle
being the segment FT of the small circle FI.
|
[top]
SMALL CIRCLE METHOD EVALUATION
The small circle subdivision method can be evaluated
in many ways, such as to what degree it satisfies the Goodchild
criteria for global grids (Kimerling et al. 1999). We have chosen,
however, to focus on two primary Goodchild criteria, area variation
and high compactness constancy. Area and compactness metrics were
computed for several competing methods, but we focus on comparisons
with the ISEA grid which has received serious consideration due
to its equal area nature and compactness consistency.
[top]
Spherical Area Variation
The area computations central to small circle
partitioning should insure that the sub-triangles will be identical
in spherical surface area at all levels of recursion. To prove that
the equation coding was correct, computed triangle areas were stored
up to 9-fold partitioning, recursion level 8 (43,046,721 triangles,
each 0.59 km2 on the globe). Statistical analysis of
these data showed a range, and hence standard deviation, of 0.0
for each level of recursion, the same result obtained for ISEA triangles.
For each level of recursion, there is no variation in triangle areas
on the sphere.
[top]
Compactness Variation
The compactness of grid cells is an important
geometrical property, as the ideal global grid would have identically
shaped cells (zero compactness variation) at the highest possible
compactness. The first way to study compactness variation is through
visual examination of grid cell appearance across the icosahedron
face. Comparing the 4-fold subdivision, recursion level 4 triangular
grid cells for the ISEA method and from the small circle subdivision
method (Figure 18), the smoother, more regular cell edges
resulting from small circle subdivision are
|
Figure 18. Small circle and Snyder
(ISEA) subdivision of a spherical icosahedron face using 4-fold
subdivision at recursion level 4. The icosahedron face is
shown on an orthographic map projection to show the earth's
sphericity.
|
readily apparent. Cell edge variations in the
ISEA grid appear to be localized along three lines radiating from
the center of the icosahedron face to each corner, bisecting the
vertex angles. Triangle edges along these lines appear to alternately
bow inward and outward relative to geodesic lines connecting two
corners. Triangle edges away from these three lines quickly smoothen
and triangles appear more equilateral in shape and similar in compactness.
Visual observation of shape and compactness variation
can be augmented by quantitative analysis of compactness variation
among cells if a suitable compactness metric exists. The Zone Standardized
Compactness (ZSC) measure, described in Kimerling et al. (1999)
was used so that comparisons could be drawn with ZSC values for
ISEA grid triangles at different recursion levels. ZSC values are
scaled relative to a 1.0 value for a spherical zone (cap) of the
same surface area as the sub-triangle.
ZSC values for 4-fold small circle subdivision
at recursion levels 1-6 were computed, with the values for recursion
levels 2,4 and 6 plotted on the equilateral triangle diagrams (not
the actual sub-triangle shapes) in Figure 19, along with corresponding
ZSC values for the ISEA method. The spatial pattern of ZSC compactness
variation for small circle subdivision at first appears chaotic,
but subtle symmetries along three lines bisecting the vertex angles
soon become apparent. At all levels of recursion the three corner
triangles, along with several interior triangles along the bisector
lines, have the highest compactness. Note that the minimum and maximum
compactness values seen on the diagrams and in the legend define
the range of compactness for the ISEA grid, but for both methods
compactness values are well below and above the 0.7776 compactness
of a spherical equilateral triangle. This is because for both methods
the triangle edges along the bisector lines bow outward and inward,
alternately increasing and decreasing triangle compactness. However,
the amount of bowing, and hence the compactness range, is much greater
with the ISEA method. The triangular diagrams for small circle subdivision
recursion levels 4 and 6 show that this
|
Figure 19. Triangular diagrams displaying
Zone Standardized Compactness values for 4-fold Small Circle
and Snyder (ISEA) subdivision of an icosahedron face at recursion
levels 2,4 and 6.
|
alternating inward and outward bowing of small
circle edges to achieve equal area sub-triangles continues at higher
recursion levels. The general pattern is that the higher compactness
sub-triangles at the next lower recursion level are partitioned
similarly, with center triangles of lower compactness surrounded
by three higher compactness outer triangles. However, the pattern
is reversed for the center triangles of lower compactness. Here
a center triangle of higher compactness is created, surrounded by
three lower compactness outer sub-triangles. This roughly fractal-like
"triangular checkerboard" pattern of sub-triangle compactness
continues at higher recursion levels, but the compactness differences
among adjacent sub-triangles lessen except at the initial icosahedron
corners where both the highest and lowest compactness sub-triangles
are always found.
Statistical analysis of compactness in spherical
triangle grid cells gives us additional ways to compare small circle
subdivision with other partitioning methods such as the ISEA. The
range in compactness for the small circle method (Figure 20)
shows a slight increase for 4-fold subdivision with increasing levels
of recursion because the inward or outward bowing of sub-triangle
edges becomes greater at higher recursion levels. However, these
ranges are less than those for the ISEA method.
|
Figure 20.
Range of zone standardized sub-triangle compactness for the
Small Circle and Snyder (ISEA) methods for different 4-fold
subdivision levels of recursion.
|
The mean ZSC at different recursion levels behaves
similarly for both 4- and 9-fold subdivision, with 4-fold subdivision
results shown in (Figure 21). There is with a rapid decrease
in ZSC to a stable minimum value at higher recursion levels. The
small circle means are higher than those for the ISEA method, but
less than the means for equilateral triangles. Similar to the ZSC
mean values, the standard deviations show a decrease for both 4-
and 9-fold subdivision to a stable minimum at higher recursion levels
(Figure 22). The minimum standard deviation reached for small
circle subdivision is lower than that for the ISEA method.
These results for the full population of triangular
cells at recursion levels 2-8 show the small circle method to perform
better both in terms of the average compactness produced and the
variation in compactness among sub-triangles.
|
Figure 21. Standard deviation
of zone standardized sub-triangle compactness for the Small
Circle and Snyder (ISEA) methods for different 4-fold subdivision
levels of recursion.
|
|
Figure 22. Mean of zone
standardized sub-triangle compactness for the Small Circle
and Snyder (ISEA) methods for different 4-fold subdivision
levels of recursion.
|
[top]
CONCLUDING REMARKS
We have conducted a complete analysis
of the small circle subdivision method, including its mathematical
derivation, numerical implementation and detailed analysis of results.
We have also performed a comprehensive comparison of the method
with five other global partitioning methods, although only the ISEA
compactness results are presented here. The numerical comparison
methods discussed here are closely linked to survey sampling requirements
for a global grid, as expressed in the Goodchild criteria (Kimerling
et al. 1999). We focused on the comparison of surface area and compactness
variations when identical subdivision frequencies and recursion
levels were applied to these partitioning methods.
The comparison results show that
although the small circle subdivision, the ISEA, and the Tobler-Chen
methods satisfy the equal-area grid cell criterion, the small circle
subdivision method creates at each level of recursion triangular
cells with a higher mean compactness, lower compactness range and
smaller standard deviation in compactness. Thus, if the equal-area
sampling units are a prime requirement, the small circle subdivision
method is the best option among the partitioning methods examined
here.
Although the small circle subdivision
method has better compactness characteristics than the ISEA or Tobler-Chen
method, its mean compactness is lower and the range and standard
deviation in compactness are large when compared to other non equal-area
partitioning methods not examined in this article, such as the Fuller-Gray
and Gnomonic map projection surfaces (White et al. 1998). This compactness
variation is more obvious at higher 9-fold subdivision recursion
levels. One of the factors contributing to the compactness variation
is the location of the inner point in 9-fold subdivision. If we
locate the inner point in 9-fold subdivision by maximizing the mean
compactness of the surrounding six sub-triangles, the compactness
characteristics could be improved, but it is computationally complex.
Computational complexity, not geometrical
regularity, appears to be the major drawback to widespread use of
the small circle subdivision method. C-language code for its implementation,
available from the authors, bears this out. However, run times for
higher recursion levels were not prohibitively longer than the ISEA
method, so perhaps the ever increasing CPU power will soon make
the small circle method computationally practical as well as geometrically
superior to other global grid alternatives. Computation speed may
also be increased significantly by converting the equations based
on spherical trigonometry into their vector algebra equivalents.
This would very likely result in a simpler code used to generate
the numerical solution.
[top]
REFERENCES
Baumgardner, J.R. and P.O. Frederickson. 1985. Icosahedral
Discretization of the Two-Sphere. S.I.A.M. Journal of Numerical
Analysis. 22(6): 1107-15.
Brooks, D.R. 1981. Grid Systems for Earth Radiation Budget Experiment
Applications. NASA Technical Memorandum 83233. 40 pp.
Coexeter, H.S.M. 1962. The Problem of Packing a Number of Equal
Nonoverlapping Circles on a Sphere. Transactions, New York Academy
of Science, Division of Mathematics, Series II, vol. 24: 320-31.
Dutton, G.H. 1999. A Hierarchical Coordinate System for Geoprocessing
and Cartography. Berlin: Springer-Verlag. 230 pp.
Goodchild, M.F. 1994. Geographical Grid Models for Environmental
Monitoring and Analysis across the Globe (panel session). GIS/LIS
'94 Conference, Phoenix, AZ.
Kahn, R. and A. Braverman 1999. What Shall We Do with the Data
we are Expecting From Upcoming Earth Observation Satellites? Journal
of Computational and Graphical Statistics, 8(3): 575-588.
Kimerling, A.J., Sahr, K., White, D. and L. Song. 1999. Comparing
Geometrical Properties of Global Grids. Cartography and Geographic
Information Science. 26(4): 271-287.
Rakhmanov, E.A., Saff, E.B., and Y.M. Zhou. 1994. Minimal Discrete
Energy on the Sphere. Mathematical Research Letters. 1:647-62.
Saff, E.B. and A.B.J. Kuijlaars. 1997. Distributing Many Points
on a Sphere. The Mathematical Intelligencer. 19(1): 5-11.
Snyder, J.P. 1992. An Equal-Area Map Projection for Polyhedral
Globes. Cartographica. 29(1):10-21.
Tobler, W.R. and Z. Chen. 1986. A Quadtree for Global Information
Storage. Geographical Analysis. 18(4): 360-71.
White, D., Kimerling, A.J., Sahr, K. and L. Song. 1998. Comparing
Area and Shape Distortion on Polyhedral-based Recursive Partitions
of the Sphere. International Journal of Geographical Information
Science. 12: 805-27.
[top]
|
