Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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

share|improve this question

1 Answer 1

up vote 6 down vote accepted

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 indicatrix

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)
  length.meridian <- radius.meridian                                  # (4-19)
  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, 
           meridian_radius=radius.meridian, meridian_length=length.meridian,
           normal_radius=radius.normal, normal_length=length.normal,
           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
       proj.in=CRS("+init=epsg:4267"),   # NAD 27
       proj.out=CRS("+init=esri:54030")) # World Robinson projection

i <- indicatrix(r, scale=10^4, n=71)
plot(i$outline, type="n", asp=1, xlab="Easting", ylab="Northing")
polygon(i$base, col=rgb(0, 0, 0, .025), border="Gray")
lines(i$d.lambda, lwd=2, col="Gray", lty=2)
lines(i$d.phi, lwd=2, col=rgb(.25, .7, .25), lty=2)
lines(i$axis.major, lwd=2, col=rgb(.25, .25, .7))
lines(i$axis.minor, lwd=2, col=rgb(.7, .25, .25))
lines(i$outline, asp=1, lwd=2)
share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.