An example Tissot ellipse for an Equirectangular Projection?

I am trying to calculate distortion of an equirectangular projection via Tissot indicatrices. I have tried following the directions on this post (among other things) but it is beyond my understanding as an amateur.

So, I am wondering if someone would be so kind as to calculate a single Tissot ellipse for an example equirectangular lat/long (whichever is your favorite and is distorted on an equirectangular projection). I don't understand exactly what the variables are and where they come from, so seeing the equations in action would be very useful.

I am trying to understand these equations so I can plug them into a mapping program I am coding. I asked a bunch of general questions in this thread, but I think a specific example will help me figure the rest out.

Many thanks, as always.

NCashew

For the record, here is a complete, commented implementation of the Tissot indicatrix (and related) calculations in `R`, with a worked example. The source of the equations is John Snyder's Map Projections--A Working Manual.

``````tissot <- function(lambda, phi, prj=function(z) z+0, asDegrees=TRUE, A = 6378137, f.inv=298.257223563, ...) {
#
# Compute properties of scale distortion and Tissot's indicatrix at location `x` = c(`lambda`, `phi`)
# using `prj` as the projection.  `A` is the ellipsoidal semi-major axis (in meters) and `f.inv` is
# the inverse flattening.  The projection must return a vector (x, y) when given a vector (lambda, phi).
# (Not vectorized.)  Optional arguments `...` are passed to `prj`.
#
# Source: Snyder pp 20-26 (WGS 84 defaults for the ellipsoidal parameters).
# All input and output angles are in degrees.
#
to.degrees <- function(x) x * 180 / pi
to.radians <- function(x) x * pi / 180
clamp <- function(x) min(max(x, -1), 1)                             # Avoids invalid args to asin
norm <- function(x) sqrt(sum(x*x))
#
# Precomputation.
#
if (f.inv==0) {                                                     # Use f.inv==0 to indicate a sphere
e2 <- 0
} else {
e2 <- (2 - 1/f.inv) / f.inv                                       # Squared eccentricity
}
if (asDegrees) phi.r <- to.radians(phi) else phi.r <- phi
cos.phi <- cos(phi.r)                                               # Convenience term
e2.sinphi <- 1 - e2 * sin(phi.r)^2                                  # Convenience term
e2.sinphi2 <- sqrt(e2.sinphi)                                       # Convenience term
if (asDegrees) units <- 180 / pi else units <- 1                    # Angle measurement units per radian
#
# Lengths (the metric).
#
radius.meridian <- A * (1 - e2) / e2.sinphi2^3                      # (4-18)
radius.normal <- A / e2.sinphi2                                     # (4-20)
length.normal <- radius.normal * cos.phi                            # (4-21)
#
# The projection and its first derivatives, normalized to unit lengths.
#
x <- c(lambda, phi)
d <- numericDeriv(quote(prj(x, ...)), theta="x")
z <- d[1:2]                                                         # Projected coordinates
names(z) <- c("x", "y")
g <- attr(d, "gradient")                                            # First derivative (matrix)
g <- g %*% diag(units / c(length.normal, length.meridian))          # Unit derivatives
dimnames(g) <- list(c("x", "y"), c("lambda", "phi"))
g.det <- det(g)                                                     # Equivalent to (4-15)
#
# Computation.
#
h <- norm(g[, "phi"])                                               # (4-27)
k <- norm(g[, "lambda"])                                            # (4-28)
a.p <- sqrt(max(0, h^2 + k^2 + 2 * g.det))                          # (4-12) (intermediate)
b.p <- sqrt(max(0, h^2 + k^2 - 2 * g.det))                          # (4-13) (intermediate)
a <- (a.p + b.p)/2                                                  # (4-12a)
b <- (a.p - b.p)/2                                                  # (4-13a)
omega <- 2 * asin(clamp(b.p / a.p))                                 # (4-1a)
theta.p <- asin(clamp(g.det / (h * k)))                             # (4-14)
conv <- (atan2(g["y", "phi"], g["x","phi"]) + pi / 2) %% (2 * pi) - pi # Middle of p. 21
#
# The indicatrix itself.
# `svd` essentially redoes the preceding computation of `h`, `k`, and `theta.p`.
#
m <- svd(g)
axes <- zapsmall(diag(m\$d) %*% apply(m\$v, 1, function(x) x / norm(x)))
dimnames(axes) <- list(c("major", "minor"), NULL)

return(list(location=c(lambda, phi), projected=z,
scale.meridian=h, scale.parallel=k, scale.area=g.det, max.scale=a, min.scale=b,
to.degrees(zapsmall(c(angle_deformation=omega, convergence=conv, intersection_angle=theta.p))),
axes=axes, derivatives=g))
}
indicatrix <- function(x, scale=1, ...) {
# Reprocesses the output of `tissot` into convenient geometrical data.
o <- x\$projected
base <- ellipse(o, matrix(c(1,0,0,1), 2), scale=scale, ...)             # A reference circle
outline <- ellipse(o, x\$axes, scale=scale, ...)
axis.major <- rbind(o + scale * x\$axes[1, ], o - scale * x\$axes[1, ])
axis.minor <- rbind(o + scale * x\$axes[2, ], o - scale * x\$axes[2, ])
d.lambda <- rbind(o + scale * x\$derivatives[, "lambda"], o - scale * x\$derivatives[, "lambda"])
d.phi <- rbind(o + scale * x\$derivatives[, "phi"], o - scale * x\$derivatives[, "phi"])
return(list(center=x\$projected, base=base, outline=outline,
axis.major=axis.major, axis.minor=axis.minor,
d.lambda=d.lambda, d.phi=d.phi))
}
ellipse <- function(center, axes, scale=1, n=36, from=0, to=2*pi) {
# Vector representation of an ellipse at `center` with axes in the *rows* of `axes`.
# Returns an `n` by 2 array of points, one per row.
theta <- seq(from=from, to=to, length.out=n)
t((scale * t(axes))  %*% rbind(cos(theta), sin(theta)) + center)
}
#
# Example: analyzing a GDAL reprojection.
#
library(rgdal)

prj <- function(z, proj.in, proj.out) {
z.pt <- SpatialPoints(coords=matrix(z, ncol=2), proj4string=proj.in)
w.pt <- spTransform(z.pt, CRS=proj.out)
return(w.pt@coords[1, ])
}
r <- tissot(130, 54, prj,                # Longitude, latitude, and reprojection function