## Documentation Center |

On this page… |
---|

Because no projection can preserve all directional and nondirectional geographic characteristics, it is useful to be able to estimate the degree of error in direction, area, and scale for a particular projection type and parameters used. Several Mapping Toolbox™ functions display projection distortions, and one computes distortion metrics for specified locations.

A standard method of visualizing the distortions introduced
by the map projection is to display small circles at regular intervals
across the globe. After projection, the small circles appear as ellipses
of various sizes, elongations, and orientations. The sizes and shapes
of the ellipses reflect the projection distortions. Conformal projections
have circular ellipses, while equal-area projections have ellipses
of the same area. This method was invented by Nicolas Tissot in the
19th century, and the ellipses are called *Tissot indicatrices* in
his honor. The measure is a tensor function of location that varies
from place to place, and reflects the fact that, unless a map is conformal,
map scale is different in every direction at a location.

As the following example illustrates, you can add the indicatrices
to a map display with the command `tissot` and remove
them with `clmo tissot`:

Set up a Sinusoidal projection in a skewed aspect, plotting the graticule:

figure; axesm sinusoid gridm on;framem on; setm(gca,'Origin', [20 30 45])

Load the

`coast`data set and plot it as green patches:load coast patchm(lat,long,'g')

Plot the default Tissot diagram, shown below:

tissot

Notice that the circles vary considerably in shape. This indicates that the Sinusoidal projection is not conformal. Despite the distortions, however, the circles all cover equal amounts of area on the map, because the projection has the equal-area property.

Default Tissot diagrams are drawn with blue unfilled 100-point circles spaced 30 degrees apart in both directions. The default circle radius is 1/10 of the current radius of the reference ellipsoid (by default that radius is

`1`).Now clear the Tissot diagram, rotate the projection to a polar aspect, and plot a new Tissot diagram using circles paced 20 degrees apart, half as big as before, drawn with 20 points, and drawn in red:

clmo tissot setm(gca,'Origin', [90 0 45]) tissot([20 20 .05 20],'Color','r')

The result is shown below. Note that circles are drawn faster because fewer points are computed for each one. Also note that the distortions are still smallest close to the map origin, and still greatest near the map frame.

Try changing the map projection to a conformal one such as Mercator or Stereographic to see what Tissot indicatrices look like on shape-preserving maps.

For further information, see the reference page for `tissot`.

Most map projection distortions are rather orderly and vary
continuously, making them suitable for display via isolines (contour
lines). In addition to Tissot diagrams, the toolbox can plot isolines
of variations of several parameters associated with map projections,
using `mdistort`.

The `mdistort` function can plot variations
in angles, areas, maximum and minimum scale, and scale along parallels
and meridians, in units of percent deviation (except for angles, for
which degrees are used). Use this function in selecting projections
and projection parameters when you are concerned about keeping specific
types of distortion within limits. Below are some examples of `mdistort` using
the Hammer modified azimuthal projections and the Bonne pseudoconic
projection.

Create a Hammer projection map axes in normal aspect, and plot a graticule, frame, and coastlines on it:

figure; axesm('MapProjection','hammer','Grid','on','Frame','on')

Load the

`coast`data set and plot it as green patches:load coast patchm(lat,long,'g')

Call

`mdistort`to plot contours of minimum-to-maximum scale ratios:mdistort('scaleratio')

Notice that the region of minimum distortion is centered around (

`0,0`).Repeat this diagram with a Bonne projection in a new figure window:

figure; axesm('MapProjection','bonne','Grid','on','Frame','on') patchm(lat,long,'g') mdistort('scaleratio')

Notice that the region of minimum distortion is centered around (

`30,0`), which is where the single standard parallel is.You can toggle the isolines by typing

`mdistort`or`mdistort off`. Look at some other types of distortion. The types you can request are

`area`— Percent departures from equal area`angles`— Angular distortion of right angles`scale`or`maxscale`— Percent of maximum scale`minscale`— Percent of minimum scale`parscale`— Percent of scale along the parallels`merscale`— Percent of scale along the meridians`scaleratio`— Percent of maximum-to-minimum scale ratio

For further information see the reference page for `mdistort`.

The `tissot` and `mdistort` functions
described above provide synoptic visual overviews of different forms
of map projection error. Sometimes, however, you need numerical estimates
of error at specific locations in order to quantify or correct for
map distortions. This is useful, for example, if you are sampling
environmental data on a uniform basis across a map, and want to know
precisely how much area is associated with each sample point, a statistic
that will vary by location and be projection dependent. Once you have
this information, you can adjust environmental density and other statistics
you collect for areal variations induced by the map projection.

A Mapping Toolbox function returns location-specific map
error statistics from the current projection or an mstruct. The `distortcalc` function
computes the same distortion statistics as `mdistort` does,
but for specified locations provided as arguments. You provide the
latitude-longitude locations one at a time or in vectors. The general
form is

[areascale,angdef,maxscale,minscale,merscale,parscale] = ... distortcalc(mstruct,lat,long)

However, if you are evaluating the current map figure, omit the mstruct. You need not specify any return values following the last one of interest to you.

The following exercise uses `distortcalc` to
compute the maximum area distortion for a map of Argentina from the
landareas data set.

Read the North and South America polygon:

Americas = shaperead('landareas','UseGeoCoords',true, ... 'Selector', {@(name) ... strcmpi(name,{'north and south america'}),'Name'});

Set the spatial extent (map limits) to contain the southern part of South America and also include an area closer to the South Pole:

mlatlim = [-72.0 -20.0]; mlonlim = [-75.0 -50.0]; [alat, alon] = maptriml([Americas.Lat], ... [Americas.Lon], mlatlim, mlonlim);

Create a Mercator cylindrical conformal projection using these limits, specify a five-degree graticule, and then plot the outline for reference:

figure; axesm('MapProjection','mercator','grid','on', ... 'MapLatLimit',mlatlim,'MapLonLimit',mlonlim,... 'MLineLocation',5, 'PLineLocation',5) plotm(alat,alon,'b')

The map looks like this:

Sample every tenth point of the patch outline for analysis:

alats = alat(1:10:numel(alat)); alons = alon(1:10:numel(alat));

Compute the area distortions (the first value returned by distortcalc) at the sample points:

adistort = distortcalc(alats, alons);

Find the range of area distortion across Argentina (percent of a unit area on, in this case, the equator):

adistortmm = [min(adistort) max(adistort)] adistortmm = 1.1790 2.7716

As Argentina occupies mid southern latitudes, its area on a Mercator map is overstated, and the errors vary noticeably from north to south.

Remove any

`NaN`s from the coordinate arrays and plot symbols to represent the relative distortions as proportional circles, using`scatterm`:nanIndex = isnan(adistort); alats(nanIndex) = []; alons(nanIndex) = []; adistort(nanIndex) = []; scatterm(alats,alons,20*adistort,'red','filled')

The resulting map is shown below:

The degree of area overstatement would be considerably larger if it extended farther toward the pole. To see how much larger, get the area distortion for 50ºS, 60ºS, and 70ºS:

a=distortcalc(-50,-60) a = 2.4203 a=distortcalc(-60,-60) a = 4 >> a=distortcalc(-70,-60) a = 8.5485

Using this technique, you can write a simple script that lets you query a map repeatedly to determine distortion at any desired location. You can select locations with the graphic cursor using

`inputm`. For example,[plat plon] = inputm(1) plat = -62.225 plon = -72.301 >> a=distortcalc(plat,plon) a = 4.6048

Naturally the answer you get will vary depending on what point you pick. Using this technique, you can write a simple script that lets you query a map repeatedly to determine any distortion statistic at any desired location.

Try changing the map projection or even the orientation vector
to see how the choice of projection affects map distortion. For further
information, see the reference page for `distortcalc`.

Was this topic helpful?