A Tissot Indicatrix is useful method for communicating at a glance the kinds of distortion a given projection is prone to (in the figure below, each of the red circles occupies the same area). I've been told that popular methods for generating TI's have their own problems, to the point of being sometimes woefully inaccurate.
What is the problem with the popular methods, and what is most correct way of generating a TI that is accessible to your average GIS dude(ette)?
Any software that can project coordinates accurately can compute accurate Tissot indicatrices.
A good source for the formulas is Snyder, John, Map Projections--A Working Manual, primarily at pp 20-26. (I won't reproduce them here because this site doesn't have appropriate tools for communicating mathematical formulas.) They require all four first derivatives of the projected coordinates (x,y) with respect to the spherical coordinates (lat, lon) = (phi, lambda):
dx / d(phi), dx / d(lambda);
dy / d(phi), dy / d(lambda).
Everything else about the TI's is computed in terms of these (using some arithmetic and trigonometric functions: the cosine, principal inverse sine, and principal inverse tangent). The computations require a description of the earth's shape. For the greatest accuracy use an ellipsoidal datum with semimajor axis a and eccentricity e. (These will be known to the software.)
Snyder's book has instructions on how to compute everything except these derivatives. Do it numerically. I have had excellent results using first-order central finite difference estimates at a distance of h = 10^(-5.2) radians (typically around 50 meters): this is a good compromise between trying to get infinitesimally close and losing too much precision from floating point roundoff (assuming double precision), because the error made is proportional to (10^(-5.2))^2 = 10^(-10.4) and 10^(-5.2) equals 10^10.4 times the IEEE double precision accuracy of 10^(-15.6) and it's still a lot larger than typical precision in projections, which usually run from 10^(-10) to about 10^(-14).
So, how do you compute finite difference estimates? This part is surprisingly easy. To obtain dx/d(phi) at a point (phi, lambda), ask your GIS to project the points
(phi - h/2, lambda) --> (x0,y0),
(phi + h/2, lambda) --> (x1,y1).
Use the estimates
dx / d(phi) = (x1 - x0)/h,
dy / d(phi) = (y1 - y0)/h.
Similarly, project the points
(phi, lambda - h/2) --> (x2,y2),
(phi, lambda + h/2) --> (x3,y3)
and use the estimates
dx / d(lambda) = (x3 - x2)/h,
dy / d(lambda) = (y3 - y2)/h.
That takes four projections and a tiny bit of arithmetic. (You can reduce it to three by using non-central differences, but the accuracy goes down a little. It's wise aim for high accuracy, without letting h get too small, unless you are sure your GIS is using survey-grade (millimeter) accuracy in its projection formulas.)
From these derivatives, along with Snyder's formulas (paying attention to the modifications described at 4-19 and 4-21), you can obtain the lengths of the axes of the Tissot Indicatrix at (phi, lambda) and its orientation. On world-scale maps the TI will be so small as to be invisible, so the final thing to do is decide how much you want to rescale each TI. I determine the scale factor by finding out how large the map will be, finding the sizes of typical TIs across the map, and scaling so that those TIs will be approximately 6% as wide as the map. It's a good start, anyway; I let the user adjust the size of the TI from there. Of course you will rescale all the TIs by the same amount, so they can be compared, and each will be rescaled around its own center (which is obtained by a fifth projection, (phi, lambda) --> (x,y)).
A nice addition to the elliptical portrayal of the TI is to show the directions of the local meridian and parallel: then, at a glance, you can assess the grid convergence. I also show a standard circle (representing no distortion) concentric with each TI because it improves the reader's ability to gauge the amount of distortion represented by each ellipse.
Of note in this Mollweide projection is the extreme TI near the south pole. It is still a perfect ellipse and accurately describes the map distortion there.
