Quantcast

Documentation Center

  • Trial Software
  • Product Updates

Visualizing and Quantifying Projection Distortions

Displays of Spatial Error in Maps

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.

Visualizing Projection Distortions via Tissot Indicatrices

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

  1. 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])
  2. Load the coast data set and plot it as green patches:

    load coast
    patchm(lat,long,'g')
  3. 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).

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

Visualizing Projection Distortions via Isolines

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.

  1. 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')
  2. Load the coast data set and plot it as green patches:

    load coast
    patchm(lat,long,'g')
  3. 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).

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

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

Quantifying Map Distortions at Point Locations

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.

Using distortcalc to Determine Map Projection Geometric Distortions

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

  1. Read the North and South America polygon:

    Americas = shaperead('landareas','UseGeoCoords',true, ...
        'Selector', {@(name) ...
        strcmpi(name,{'north and south america'}),'Name'});
  2. 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);
  3. 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:

  4. Sample every tenth point of the patch outline for analysis:

    alats = alat(1:10:numel(alat));
    alons = alon(1:10:numel(alat));
  5. Compute the area distortions (the first value returned by distortcalc) at the sample points:

    adistort = distortcalc(alats, alons);
  6. 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.

  7. Remove any NaNs 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:

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

  9. 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?