In mdsumner/tissot: The Tissot Indicatrix

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.path = "man/readmefigs/README-", 
        fig.height = 7, 
       fig.width = 7
)

The Tissot Indicatrix

The Tissot Indicatrix is used to characterize local distortions within map projections.

I have derived the code in this package (with permission) from Bill Huber's wonderful online answer here:

http://gis.stackexchange.com/questions/31651/an-example-tissot-ellipse-for-an-equirectangular-projection

Also see

https://gis.stackexchange.com/questions/5068/how-to-create-an-accurate-tissot-indicatrix

Installation

Can be installed with

remotes::install_github("hypertidy/tissot")

Minimal example

library(tissot)
# NAD 27 in
# World Robinson projection out
r <- tissot(130, 54,
           proj.in= "EPSG:4267",  
           proj.out= "ESRI:54030")
i0 <- indicatrix(r, scale=1e2, n=71)
plot(i0)
plot(tissot:::.prj(earthcircle::earthcircle(130, 54, scale = 5e2), tissot_get_proj()))
tissot_abline(130, 54)

tissot_map(add = FALSE, xlim = c(8.5e6, 1.3e7), ylim = c(4e6, 7e6))
i1 <- indicatrix(r, scale=1e6, n=71)
plot(i1, add = T)
tissot_abline(130, 54)
lines(tissot:::.prj(earthcircle::earthcircle(130, 54, scale = 2e6), tissot_get_proj()) %*% (diag(2) * 2))

Since an original port of whuber's code we have now made it much easier to create many indicatrixes and plot them in one step. Or we can still just grab one and plot it on its own. Note that the scale is quite different in these plots.

x <- seq(-172.5, 172.5, by = 15)
y <- seq(-82.5, 82.5, by = 15)
xy <- expand.grid(x, y)
r <- tissot(xy,
            proj.in= "OGC:CRS84",
            proj.out= "+proj=robin")

j <- which.min(abs(135 - r$lon) + abs(54 - r$lat))
i <- indicatrix0(r[j, ], scale= 1e4, n=71)
plot(i, add = FALSE)

ii <- indicatrix(r, scale=4e5, n=71)
tissot_map(add = FALSE)
plot(ii, add = TRUE)
tissot_abline(xy[j, 1], xy[j, 2])

Mollweide.

m <- tissot(xy,
            proj.in= "OGC:CRS84",
            proj.out= "+proj=moll")


plot(indicatrix(m, scale=4e5, n=71), add = FALSE)
tissot_map()

Eckhert III

e <- tissot(xy,
            proj.in= "OGC:CRS84",
            proj.out= "+proj=eck3")


plot(indicatrix(e, scale=4e5, n=71), add = FALSE)

Equidistant

aeqd <- tissot(xy,
            proj.in= "OGC:CRS84",
            proj.out= "+proj=aeqd")


plot(indicatrix(aeqd, scale=4e5, n=71), add = FALSE)

Cassini-Soldner (spherical because ellipsoidal seems broken)

xx <- tissot(xy,
            proj.in= "OGC:CRS84",
            proj.out= "+proj=cass +R=6378137")


plot(indicatrix(xx, scale=4e5, n=71), add = FALSE)
points(tissot_map(col = "transparent"), pch = ".")

Sinusoidal

s <- tissot(xy,
            proj.in= "OGC:CRS84",
            proj.out= "+proj=sinu")


plot(indicatrix(s, scale=3e5, n=71), add = FALSE)
tissot_map()

Polar example

p <- tissot(xy[xy[,2] < -30, ],
            proj.in= "OGC:CRS84",
            proj.out= "+proj=stere +lon_0=147 +lat_ts-71 +lat_0=-90 +datum=WGS84")

plot(indicatrix(p, scale = 3e5))
tissot_map()
tissot_abline(147, -42)

laea <- tissot(xy[xy[,2] < 20, ],
            proj.in= "OGC:CRS84",
            proj.out= "+proj=laea +lon_0=147 +lat_0=-90 +datum=WGS84")

plot(indicatrix(laea, scale = 3e5))

Oblique Mercator

You get the idea ... many projections need extra attention for real data.

mp0 <- do.call(cbind, maps::map(plot = FALSE)[1:2])
omerc <- "+proj=omerc +lonc=147 +gamma=9 +alpha=9 +lat_0=-80 +ellps=WGS84"
mp <- tissot:::.prj(mp0, omerc, proj.in = "OGC:CRS84")
o <- tissot(xy,
            proj.in= "OGC:CRS84",
            proj.out= omerc)

plot(indicatrix(o, scale = 3e5))
lines(mp)

VicGrid

vgrid <- "+proj=lcc +lat_1=-36 +lat_2=-38 +lat_0=-37 +lon_0=145 +x_0=2500000 +y_0=2500000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
mp <- tissot:::.prj(mp0, vgrid, proj.in = "OGC:CRS84")
v <- tissot(as.matrix(expand.grid(seq(120, 165, by =5 ), 
                                          seq(-45, -35, by = 5))),
            proj.in= "OGC:CRS84",
            proj.out= vgrid)

plot(indicatrix(v, scale = 2e5))
lines(mp)

UTM Zone 54 (Hobart)

utm <- "+proj=utm +zone=54 +south"
mp <- tissot:::.prj(mp0, utm, proj.in = "OGC:CRS84")
u <- tissot(as.matrix(expand.grid(seq(108, 162, by =6 ), 
                                          seq(-65, 55, by = 15))),
            proj.in= "OGC:CRS84",
            proj.out= utm)

plot(indicatrix(u, scale = 2e5))
lines(mp)

library(tissot)
library(maptools)
library(raster)
buildandplot <- function(data, scale = 5e5, ...) {
  ## grid of points
  gr <- rasterToPoints(raster(data, nrow = 7, ncol = 7), spatial = FALSE)
  ## relying on dev {PROJ} that links to unreleased {libproj}
  grll <- reproj::reproj_xy(gr, "OGC:CRS84", source = projection(data) ))
  sp::plot(data,  ...)
  grll <- grll[!is.na(grll[,1]), ]
  for (i in seq_len(nrow(grll))) {
    tis <- tissot(grll[i, 1], grll[i, 2],  
                                               proj.in = projection(wrld_simpl), proj.out = projection(data))
   ind <- indicatrix(tis, scale = scale, n = 71)
   plot(ind, add = TRUE)
  }
  invisible(NULL)
}


## choose a projection
ptarget1 <- "+proj=stere +lon_0=147 +lat_ts-71 +lat_0=-90 +ellps=WGS84"
w1 <- spTransform(subset(wrld_simpl, coordinates(wrld_simpl)[,2] < 10), CRS(ptarget1))

ptarget2 <- "+proj=laea +lon_0=147 +lat_0=-90 +ellps=WGS84"
w2 <- spTransform(subset(wrld_simpl, coordinates(wrld_simpl)[,2] < 10), CRS(ptarget2))

ptarget3 <- "+proj=omerc +lonc=147 +gamma=9 +alpha=9 +lat_0=-80 +ellps=WGS84"
w3 <- spTransform(subset(wrld_simpl, coordinates(wrld_simpl)[,2] < -12), CRS(ptarget3), scale = 3e5)


buildandplot(w1, main = "Polar Stereographic")
buildandplot(w2, main = "Lambert Azimuthal Equal Area")

buildandplot(w3, main = "Oblique Mercator")

## doesn't look right
# ptarget8 <- "+proj=laea +lat_0=-90"
# w8 <- spTransform(wrld_simpl, CRS(ptarget8))
# buildandplot(w8)




library(raster)
ptarget4 <- "+proj=merc +ellps=WGS84"
w4 <- spTransform(raster::intersect(disaggregate(wrld_simpl), as(extent(-180, 180, -85, 90), "SpatialPolygons")), ptarget4)
buildandplot(w4, main = "Mercator")

ptarget5 <- "+proj=lcc +ellps=WGS84 +lon_0=134 +lat_0=-30 +lat_1=-50 +lat_2=-20"
w5 <- spTransform(raster::intersect(disaggregate(wrld_simpl), as(extent(80, 180, -65, -10), "SpatialPolygons")), ptarget5)
buildandplot(w5, main = "Lambert Conformal Conic", scale = 3.5e5)


ptarget6 <- "+proj=utm +zone=50 +south +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs "

w6 <- spTransform(raster::intersect(disaggregate(wrld_simpl), as(extent(80, 160, -65, -10), "SpatialPolygons")), ptarget6)
buildandplot(w6, main = "UTM South Zone 50 ", col = "grey", scale = 2.5e5)



buildandplot(wrld_simpl, main = "Longitude / Latitude")
degAxis(1)
degAxis(2)

## changes in spatial break this old hack
library(dplyr)
ex <- extent(c(20891678,  40158321, -13438415,  10618277))
target7 <- "+proj=lsat +lsat=5 +path=188"
library(spbabel)
tab <- sptable(spTransform(disaggregate(wrld_simpl), target7)) %>% filter(x_ >= xmin(ex), x_ <= xmax(ex), y_ >= ymin(ex), y_ <= ymax(ex))
## egregiously naive crop here, but good enough for the task
w7 <- sp(tab  %>% group_by(branch_)  %>% summarize(n = n())  %>% filter(n > 2) %>% inner_join(tab), crs = target7)
library(graticule)
g <- graticule(seq(-180, 165, by = 15), seq(-85, -20, by = 5), proj = target7, xlim = c(-180, 180), ylim = c(-85, -5))
buildandplot(w7, main = "Space Oblique Mercator, lsat=5, path=188 ", col = "grey", scale = 5e5)
plot(g, add = TRUE, lty = 2)