R-CMD-check

tissot

The Tissot Indicatrix characterizes local distortion in map projections. This package computes and plots indicatrixes using a vectorized finite-difference Jacobian with gdalraster as the projection engine.

Derived (with permission) from Bill Huber’s GIS StackExchange answer.

Installation

# install.packages("pak")
pak::pak("hypertidy/tissot@refactor-2026")

Quick start

tissot() returns a tibble of distortion properties. The second argument is the projection target; source defaults to EPSG:4326:

library(tissot)
tissot(c(147, -42), "+proj=utm +zone=55 +south")
#> Tissot indicatrix: 1 point, +proj=utm +zone=55 +south
#> # A tibble: 1 × 14
#>       x     y dx_dlam    dy_dlam dx_dphi dy_dphi scale_h scale_k scale_omega
#>   <dbl> <dbl>   <dbl>      <dbl>   <dbl>   <dbl>   <dbl>   <dbl>       <dbl>
#> 1   147   -42 0.99960 -5.8386e-7       0 0.99960 0.99960 0.99960 0.000033471
#> # ℹ 5 more variables: scale_a <dbl>, scale_b <dbl>, scale_area <dbl>,
#> #   angle_deformation <dbl>, convergence <dbl>

Columns include: scale_h (meridional), scale_k (parallel), scale_a / scale_b (max/min singular values), scale_area, angle_deformation, and convergence.

Plotting indicatrixes

indicatrix() builds plottable ellipses. The dashed circle is the undistorted reference; the filled ellipse shows the projection’s distortion.

xy <- expand.grid(x = seq(0, 1e6, length.out = 5), y = seq(4900000, 5700000, length.out = 4))
lonlat <- tissot_unproject(xy, source = "+proj=utm +zone=55 +south")
tis <- tissot(lonlat, "+proj=utm +zone=55 +south")
plot(indicatrix(tis), scale = 3e4)
tissot_map()

#> GDAL FAILURE 1: Point outside of projection domain
#> GDAL FAILURE 1: Point outside of projection domain
#> GDAL FAILURE 1: Point outside of projection domain
#> GDAL FAILURE 1: Point outside of projection domain
#> GDAL FAILURE 1: Point outside of projection domain
#> GDAL FAILURE 1: Point outside of projection domain
#> GDAL FAILURE 1: Point outside of projection domain
#> GDAL FAILURE 1: Point outside of projection domain

What does that top left indicatrix look like?

plot(indicatrix(tis)[1])

Far from our UTM zone we are in a lot more trouble.

## UTM zone 55 is at 147 longitude (55 * 6 - 183)
tis <- tissot(cbind(100, -42), "+proj=utm +zone=55 +south")
plot(indicatrix(tis))


##  In Mercator we have well known problems 
tis <- tissot(cbind(147, -42), "+proj=merc")
plot(indicatrix(tis))


## close to the equator Mercator is ok (in exactly the same way that UTM Zone 55 is ok near 147E longitude)
tis <- tissot(cbind(147, 0), "+proj=merc")
plot(indicatrix(tis))

Map projection is arbitrary.

xy <- expand.grid(seq(-150, 150, by = 30), seq(-60, 60, by = 30))
r <- tissot(xy, "+proj=robin")
ii <- indicatrix(r)
plot(ii, scale = 6e5, add = FALSE, show.axes  = TRUE, show.circle = TRUE)
tissot_map()

Distortion summary

summary(r)
#> Tissot indicatrix: 55 points
#>   Source CRS: EPSG:4326
#>   Target CRS: +proj=robin
#>   Areal scale:  min=0.8209  max=1.2790  mean=1.0468
#>   Angular def:  min=1.5058  max=52.3188  mean=21.3786 deg
#>   Scale h:      min=0.8790  max=1.3023  (meridional)
#>   Scale k:      min=0.8487  max=1.3521  (parallel)

Colour-coded distortion

Pass fill.by to colour ellipses by a distortion metric:

plot(ii, scale = 6e5, add = FALSE, fill.by = "scale_area")
tissot_map()

plot(ii, scale = 6e5, add = FALSE, fill.by = "angle_deformation")
tissot_map()

Projection comparison

m <- tissot(xy, "+proj=moll")
plot(indicatrix(m), scale = 5e5, add = FALSE)
tissot_map()

merc_xy <- expand.grid(seq(-150, 150, by = 30), seq(-75, 75, by = 15))
me <- tissot(merc_xy, "+proj=merc")
plot(indicatrix(me), scale = 5e5, add = FALSE)
tissot_map()

Rich single-indicatrix plots

A single indicatrix with axes and reference circle:

ii2 <- indicatrix(c(147, -42), "+proj=lcc +lat_1=-36 +lat_2=-38 +lat_0=-37 +lon_0=145")
plot(ii2[[1]], scale = 1e4, add = FALSE, show.axes = TRUE, show.circle = TRUE)

ii3 <- indicatrix(c(147, -42), "+proj=lcc +lat_1=-15 +lat_2=5 +lat_0=0 +lon_0=145")
plot(ii3[[1]], scale = 1e4, add = FALSE, show.axes = TRUE, show.circle = TRUE)

Arbitrary projections

Lambert Conformal Conic

pxy <- expand.grid(seq(100, 200, by = 25), seq(-75, -45, by = 10))
p <- tissot(pxy, "+proj=lcc +lat_0=-60 +lon_0=147 +lat_1=-70 +lat_2=-55")
plot(indicatrix(p), scale = 3e5, add = FALSE, fill.by = "scale_area")
tissot_map()

Universal Transverse Mercator

qxy <- expand.grid(seq(100, 200, by = 25), seq(-75, -45, by = 10))
p <- tissot(qxy, "EPSG:32755")
plot(indicatrix(p), scale = 3e5, add = FALSE, fill.by = "scale_area")
tissot_map()

#> GDAL FAILURE 1: Point outside of projection domain
#> GDAL FAILURE 1: Point outside of projection domain
#> GDAL FAILURE 1: Point outside of projection domain
#> GDAL FAILURE 1: Point outside of projection domain
#> GDAL FAILURE 1: Point outside of projection domain
#> GDAL FAILURE 1: Point outside of projection domain
#> GDAL FAILURE 1: Point outside of projection domain
#> GDAL FAILURE 1: Point outside of projection domain

Polar projections

In any projection we should refer to a regular grid of points in its crs, else we get weird situations like this, more obvious on an actual pole:

polar_xy <- expand.grid(seq(-180, 150, by = 30), seq(-80, -50, by = 10))
p <- tissot(polar_xy, "+proj=stere +lat_0=-90 +lon_0=147")
plot(indicatrix(p), scale = 2.5e5, add = FALSE, fill.by = "scale_area")
tissot_map()

la <- tissot(polar_xy, "+proj=laea +lat_0=-90 +lon_0=147")
plot(indicatrix(la), scale = 2.5e5, add = FALSE)
tissot_map()

If we push away from the pole in Lambert Azimuthal Equidistant it’s useful to see what happens.

lea <- tissot(polar_xy, "+proj=aeqd +lat_0=-20 +lon_0=147")
plot(indicatrix(lea), scale = 2.5e5, add = FALSE)
tissot_map()

Consider generating input in the crs you are assessing

As with the UTM example above with tissot_unproject() it is usually far better is to generate a grid in the crs we are assessing. A grid in lon/lat won’t be very meaningful in many projections depending on the context.

op <- par(mfrow = c(1, 2))
ext <- c(-180, 150, -80, -50)
crs <- "+proj=stere +lat_0=-90 +lon_0=147"
projext <- gdalraster::bbox_transform(ext[c(1, 3, 2, 4)], srs_to = crs, srs_from = "EPSG:4326")

polar <- expand.grid(seq(projext[1L], projext[3L], by = 30 * 1e5), seq(projext[2], projext[4], by = 10 * 1e5))
polar_xy <- tissot_unproject(polar, "EPSG:4326", source = crs)
p <- tissot(polar_xy, crs, source = "EPSG:4326")
plot(indicatrix(p), scale = 2.5e5, add = FALSE, fill.by = "scale_area")
tissot_map()

ext <- c(-180, 150, -80, -50)
crs <- "+proj=laea +lat_0=-90 +lon_0=147"
projext <- gdalraster::bbox_transform(ext[c(1, 3, 2, 4)], srs_to = crs, srs_from = "EPSG:4326")

polar <- expand.grid(seq(projext[1L], projext[3L], by = 30 * 1e5), seq(projext[2], projext[4], by = 10 * 1e5))
polar_xy <- tissot_unproject(polar, "EPSG:4326", source = crs)
p <- tissot(polar_xy, crs, source = "EPSG:4326")
plot(indicatrix(p), scale = 2.5e5, add = FALSE, fill.by = "scale_area")
tissot_map()

par(op)

Why this package?

Most “Tissot indicatrix” plots you’ll find online are just geographic circles drawn on the map. They show what happens to a circle under the projection, which is useful — but it’s not the indicatrix. The indicatrix is the Jacobian of the projection at a point: it gives you actual scale factors, angular deformation, and areal distortion. This package computes those.

Other examples: mgimond.

Code of Conduct

Please note that the tissot project is released with a Contributor Code of Conduct. By contributing to this project, you agree to abide by its terms.