Documentation |
On this page… |
---|
Visualize Spatial Error in Maps Visualize Projection Distortions using Isolines |
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 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:
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.