Spheres, Spheroids, and Geoids

Working with geospatial data involves geographic concepts (e.g., geographic and plane coordinates, spherical geometry) and geodetic concepts (such as ellipsoids and datums). This group of sections explain, at a high level, some of the concepts that underlie geometric computations on spherical surfaces.

Although the Earth is very round, it is an oblate spheroid rather than a perfect sphere. This difference is so small (only one part in 300) that modeling the Earth as spherical is sufficient for making small-scale (world or continental) maps. However, making accurate maps at larger scale demands that a spheroidal model be used. Such models are essential, for example, when you are mapping high-resolution satellite or aerial imagery, or when you are working with coordinates from the Global Positioning System (GPS). This section addresses how Mapping Toolbox software accurately models the shape, or figure, of the Earth and other planets.

 Back to Top

Geoid and Ellipsoid

Literally, geoid means Earth-shaped. The geoid is an empirical approximation of the figure of the Earth (minus topographic relief), its "lumpiness.." Specifically, it is an equipotential surface with respect to gravity, more or less corresponding to mean sea level. It is approximately an oblate ellipsoid, but not exactly so because local variations in gravity create minor hills and dales (which range from -100 m to +60 m across the Earth). This variation in height is on the order of one percent of the differences between the semimajor and semiminor ellipsoid axes used to approximate the Earth's shape, as described in The Ellipsoid Vector.

Mapping the Geoid

The following figure, made using the geoid data set, maps the figure of the Earth. To execute these commands, select them all by dragging over the list in the Help browser, then click the right mouse button and choose Evaluate Selection:

load geoid; load coast

figure; axesm robinson

geoshow(geoid,geoidlegend,'DisplayType','texturemap')

colorbar('horiz')

geoshow(lat,long,'color','k')

The shape of the geoid is important for some purposes, such as calculating satellite orbits, but need not be taken into account for every mapping application. However, knowledge of the geoid is sometimes necessary, for example, when you compare elevations given as height above mean sea level to elevations derived from GPS measurements. Geoid representations are also inherent in datum definitions.

You can define ellipsoids in several ways. They are usually specified by a semimajor and a semiminor axis, but are often expressed in terms of a semimajor axis and either inverse flattening (which for the Earth, as mentioned above, is one part in 300) or eccentricity. Whichever parameters are used, as long as an axis length is included, the ellipsoid is fully constrained and the other parameters are derivable. The components of an ellipsoid are shown in the following diagram.

The toolbox includes ellipsoid models that represent the figures of the Sun, Moon, and planets, as well as a set of the most common ellipsoid models of the Earth.

 Back to Top

The Ellipsoid Vector

• Mapping Toolbox Ellipsoid Management
• Functions that Define Ellipsoid Vectors
• What Is the "Correct" Ellipsoid Vector?

Mapping Toolbox ellipsoid representations are two-element vectors, called ellipsoid vectors. The ellipsoid vector has the form [semimajor_axis eccentricity]. The semimajor axis can be in any unit of distance; the choice of units typically drives the units used for distance outputs in the toolbox functions. Meters, kilometers, or Earth radii (i.e., a unit sphere) are most frequently used. See Functions that Define Ellipsoid Vectors for details.

Eccentricity can range from 0 to 1. Most toolbox functions accept a scalar in place of an ellipsoid vector. In this case, its value is interpreted as the radius of a reference sphere, which is equivalent to an ellipsoid with an eccentricity of zero.

Standard values for the ellipsoid vector, along with several other kinds of planetary data for each of the planets and the Earth's moon, are provided by the Mapping Toolbox almanac function (see Planetary Almanac Data). In the almanac function, the default ellipsoid for the Earth is the 1980 Geodetic Reference System ellipsoid:

format long g

almanac('earth','ellipsoid','kilometers')

ans =

6378.137 0.0818191910428158

Compare this to a spherical ellipsoid definition:

almanac('earth','sphere','kilometers')

ans =

6371 0

You should set format to long g, as above, if you want to display eccentricity values at full precision.

For example, examine the parameters of the wgs72 (the 1972 World Geodetic System) ellipsoid, using the almanac function:

wgs72 = almanac('earth','wgs72','kilometers')

wgs72 =

6378.135 0.0818188106627487

Compare this with Bessel's 1841 ellipsoid:

format long g

bessel = almanac('earth','bessel','kilometers')

bessel =

6377.397155 0.0816968312225275

The ellipsoid vector's values are the semimajor axis, in kilometers, and eccentricity. Both eccentricity and flattening are dimensionless ratios. The toolbox has functions to convert elliptical definitions from these forms to ellipsoid vector form. For example, the function axes2ecc returns an eccentricity when given semimajor and semiminor axes as arguments.

The ellipse in the previous diagram is highly exaggerated. For the Earth, the semimajor axis is about 21 kilometers longer than the semiminor axis. Use the almanac function to verify this:

grs80 = almanac('earth','ellipsoid','kilometers')

grs80 =

6378.137 0.0818191910428158

semiminor = minaxis(grs80)

semiminor =

            6356.75231414036

semidiff = grs80(1) - semiminor

semidiff =

21.3846858596444

When compared to the semimajor axis, which is almost 6400 kilometers, this difference seems insignificant and can be neglected for world and other small-scale maps. For example, the scale at which 21.38 km would be smaller than a 0.5 mm line on a map (which is a typical line weight in cartography) is

kmtomm = unitsratio('mm','km')

kmtomm =

1000000

scalelim = semidiff * kmtomm / 0.5

scalelim =

4.2769e+007

The unitsratio function was used to convert the distance semidiff from kilometers into millimeters. This indicates that the Earth's eccentricity is not geometrically meaningful at scales of less than 1:43,000,000, which is roughly the scale of a world map shown on a page of this document. Consequently, most Mapping Toolbox functions default to a spherical model of the Earth. Another reason for defaulting to a sphere is that angular distances are not meaningful on ellipsoids, and some Mapping Toolbox functions compute or use angular distances. See Working with Distances on the Sphere for more information. Regardless, you are free to specify any ellipsoid when you define map axes or otherwise operate on geodata.

Mapping Toolbox Ellipsoid Management

Most maps you make with the toolbox are displayed in a map axes, which is a MATLAB axes that contains a key data structure called a "map projection structure," or mstruct. A reference ellipsoid is fundamental to defining a map axes, and is stored in the geoid field of the mstruct. (The geographic term "geoid" actually refers to a model of the shape of the earth that is much more detailed. See Geoid and Ellipsoid for more information.) Other mstruct fields specify parameters that define the map axes' current projection and for controlling the appearance of the map frame, grid, and grid labels. You define an mstruct with the axesm or defaultm functions. See Map Axes Object Properties for definitions of the fields found in mstructs.

You can pass an mstruct to certain functions you call. Other functions obtain the mstruct from the current map axes. (If it is not a map axes, such functions error.) When axesm or defaultm create a map axes containing an mstruct, their default behavior is to use a unit sphere for the ellipsoid vector. Unless you override this default, you must work in units of earth radii (or radii of whatever planet you are mapping). The following short example shows this clearly (getm obtains mstruct parameters from a map axes):

worldmap australia

ellipsoid = getm(gca,'geoid')

ans =

1 0

The worldmap function chooses map projections and parameters appropriate to the region specified to it and sets up default values for the rest of the mstruct. The geoid parameter is the ellipsoid vector that worldmap generated. The first element of the output vector indicates that the semimajor axis has a length of 1; the second element indicates that there is no eccentricity. Therefore, you are working with a sphere—a unit sphere, to be specific.

If, instead of using default ellipsoid vectors, you prefer to be explicit about your reference ellipsoid, then you can work in the length units of your choice, on either a sphere or an ellipsoid. In following example (on the sphere),

axesm('mapprojection','mercator',...

'geoid',almanac('earth','radius','meters'))

[x, y] = mfwdtran(0,90)

x =

1.0008e+07

y =

0

the projected map coordinates for a point at 0 degrees latitude, 90 degrees longitude falls just over 10⁷ meters east of the origin. If you then revert to a unit sphere (the default ellipsoid), the distance units are quite different:

axesm mercator

[x, y] = mfwdtran(0,90)

x =

1.5708

y =

0

This value for x turns out to equal π/2, which might tempt you to think that the Mercator projection has simply converted degrees to radians. But what has actually changed is that the point at (0, 90) now maps to a point 1 earth radius east of the origin. Because Mercator is a cylindrical projection having no length distortion along the equator, and because a radian is defined in terms of a sphere's radius, the numbers just happen to work out this way.

Functions that Define Ellipsoid Vectors

Some functions define a radius or an ellipsoid and can make different choices when doing so. In addition to axesm and defaultm, which create mstructs with ellipsoid vectors that default to a unit sphere, the following functions have default ellipsoid vectors or radii:

The elevation Function.  The elevation function uses the GRS 80 ellipsoid in meters as its default; unless you specify a reference ellipsoid vector yourself, elevation will assume that input altitudes and the output slant range are both in units of meters.

The distance and reckon Functions.  These functions assume by default a reference sphere with a radius of 1 (a unit sphere), but scale their range inputs and outputs to equal the size (in degrees) of the angle subtended by rays joining the center of the Earth (or planet) to the start and end points. To obtain results on an ellipsoid you must specify an ellipsoid vector such as almanac provides.

Angle-Distance Conversion Functions.  The default behavior of the 12 angle-distance conversion utilities (itemized in Working with Distances on the Sphere) is different than the above; as discussed below, these functions assume a sphere with a radius of 6371 kilometers (or, equivalently, 3440.065 nautical miles or 3958.748 statute miles), which is a reasonable average radius for Earth.

See the documentation for individual functions if you are not clear whether or how they may generate default reference ellipsoids.

What Is the "Correct" Ellipsoid Vector?

Many different reference ellipsoids have been proposed through the years. They differ because of the surveying information upon which they are based, or because they are intended to approximate the Earth only within a specific geographic region. In many cases you will want to use either the Geodetic Referencing System of 1980 (GRS80) ellipsoid or the World Geodetic System 1984 (WGS84); their semimajor axis lengths are equal and their semiminor axes (i.e., center to pole) differ in length by just over 1/10 mm, as the following code demonstrates:

grs80 = almanac('earth','grs80','meters');

wgs84 = almanac('earth','wgs84','meters');

minaxis(wgs84) - minaxis(grs80)

ans =

1.0482e-004

The toolbox supports several other ellipsoid vectors, for models ranging from Everest's 1830 ellipsoid (used for India) to the International Astronomical Union ellipsoid of 1965 (used for Australia). These can be referenced by the following statements:

ellipsoid1 = almanac('earth','ellipsoid','kilometers','everest');

ellipsoid2 = almanac('earth','ellipsoid','kilometers','iau65');

See the reference page for the almanac function for more information on the ellipsoids that are built into the toolbox. If you cannot find the ellipsoid vector you need, you can create it in the following form:

ellipsoidvec = [semimajor_axis eccentricity]