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.

This example shows how to visualize map projection distortions using isolines (contour lines). Since distortions are rather orderly and vary continuously, they are well-suited for isolines. 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).

Create a Hammer projection map axes in normal aspect and plot a graticule and frame.

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')

Plot contours of minimum-to-maximum scale ratios, using `mdistort`

. Notice that the region of minimum distortion is centered around (0,0).

```
mdistort('scaleratio')
```

Repeat this diagram with a Bonne projection in a new figure window. 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`

.

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

This example shows how to add Tissot indicatrices to a map display.

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. 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).

tissot

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. In the result, note that the 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.

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

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

**Note**You can only use`distortcalc`

to query locations that are within the current map frame or mstruct limits. Outside points yield`NaN`

as a result.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?