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)
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)
Best Answer
For this purpose simple approximations are more than good enough. North or south, one degree is about 69 miles but east or west, it is only 69*cos(latitude) miles. Because latitudes do not change much over a ten mile span, you can safely use the cosine of the central latitude of the "square." Therefore the desired coordinates for square vertices at distance r miles from a central location (f,l), given as lat-lon, are computed as
For example, suppose r = 10 miles and the central location is at latitude 50 degrees north, longitude 1 degree west, so that (f,l) = (50,-1) degrees. Then
and the coordinates are (49.855,-1.225), (50.145,-1.225), (50.145, -0.775), and (49.855, -0.775) as you march clockwise around the square starting at its southwestern corner.
Don't use this approximation near the poles or for squares larger than a few degrees on a side. Also, depending on limitations of the GIS, some care might be needed around the global cut in longitude, usually taken at +-180 degrees.