| Title: | Tissot Indicatrix for Map Projection Distortion |
|---|---|
| Description: | Compute and visualize the 'Tissot Indicatrix' for map projections. The indicatrix characterizes projection distortion by computing scale factors, angular deformation, areal distortion, and convergence at arbitrary points. Based on the calculations shared by Bill Huber on <https://gis.stackexchange.com/a/5075/482>. Uses 'PROJ' for coordinate transformation and distortion factor computation. Developed using the method published in Snyder, JP (1987) <doi:10.3133/pp1395>. |
| Authors: | Michael Sumner [aut, cre, cph] (ORCID: <https://orcid.org/0000-0002-2471-7511>), Bill Huber [aut] (original algorithm and calculations) |
| Maintainer: | Michael Sumner <[email protected]> |
| License: | MIT + file LICENSE |
| Version: | 0.3.0 |
| Built: | 2026-07-12 15:44:02 UTC |
| Source: | https://github.com/hypertidy/tissot |
Subset an indicatrix_list
## S3 method for class 'indicatrix_list' x[i]## S3 method for class 'indicatrix_list' x[i]
x |
an |
i |
indices |
An indicatrix_list with the selected elements
Computes an indicatrix_list containing per-point distortion objects
suitable for plotting with plot.indicatrix_list() or querying directly.
indicatrix(x, target = NULL, ..., source = "EPSG:4326")indicatrix(x, target = NULL, ..., source = "EPSG:4326")
x |
a |
target |
target projection CRS (extracted from |
... |
passed to |
source |
source CRS (default |
indicatrix() accepts either:
A tissot_tbl object (from tissot()) — projection is extracted
from attributes, target is optional
Any xy-ish input with an explicit target
An indicatrix_list object (a list of indicatrix objects with
source and target stored as attributes)
tissot(), plot.indicatrix(), plot.indicatrix_list(),
ti_ellipse()
## From a tissot_tbl r <- tissot(cbind(seq(-150, 150, by = 30), 0), "+proj=robin") ii <- indicatrix(r) plot(ii, scale = 5e5) ## From raw coordinates ii2 <- indicatrix(c(0, 45), "+proj=stere +lat_0=90") plot(ii2)## From a tissot_tbl r <- tissot(cbind(seq(-150, 150, by = 30), 0), "+proj=robin") ii <- indicatrix(r) plot(ii, scale = 5e5) ## From raw coordinates ii2 <- indicatrix(c(0, 45), "+proj=stere +lat_0=90") plot(ii2)
Draws a single Tissot indicatrix ellipse on the current plot. The ellipse shows the distortion of a unit circle under the map projection. Optional overlays include a reference unit circle, and lambda/phi direction axes.
## S3 method for class 'indicatrix' plot( x, scale = 1e+05, n = 72, col = "#FF990055", border = "black", add = TRUE, show.axes = TRUE, show.circle = TRUE, ... )## S3 method for class 'indicatrix' plot( x, scale = 1e+05, n = 72, col = "#FF990055", border = "black", add = TRUE, show.axes = TRUE, show.circle = TRUE, ... )
x |
an |
scale |
scaling factor for the ellipse size in projected units |
n |
number of points on the ellipse |
col |
fill colour for the ellipse |
border |
border colour |
add |
logical; add to existing plot (default |
show.axes |
|
show.circle |
|
... |
passed to |
show.circle and show.axes accept TRUE (use defaults), FALSE
(hide), or a named list of graphical parameters to override defaults.
For example, show.circle = list(border = "blue", lty = 3).
The input x, invisibly.
indicatrix(), plot.indicatrix_list(), ti_ellipse()
Draws all indicatrixes in an indicatrix_list, optionally creating a
new plot or adding to an existing one. Can colour-code the fill by a
distortion metric.
## S3 method for class 'indicatrix_list' plot( x, scale = 1e+05, n = 72, col = "#FF990055", border = "black", add = FALSE, show.axes = TRUE, show.circle = TRUE, fill.by = NULL, palette = NULL, ncolors = 64L, ... )## S3 method for class 'indicatrix_list' plot( x, scale = 1e+05, n = 72, col = "#FF990055", border = "black", add = FALSE, show.axes = TRUE, show.circle = TRUE, fill.by = NULL, palette = NULL, ncolors = 64L, ... )
x |
an |
scale |
scaling factor for ellipse size in projected units |
n |
number of points per ellipse |
col |
fill colour. If a single colour, used for all ellipses. If
|
border |
border colour |
add |
logical; add to existing plot? If |
show.axes |
|
show.circle |
|
fill.by |
character; name of a distortion metric to colour-code
the fill. One of |
palette |
colour palette function (default |
ncolors |
number of colours in the palette (default 64) |
... |
passed to |
The input x, invisibly.
indicatrix(), plot.indicatrix(), tissot_map()
xy <- expand.grid(seq(-150, 150, by = 30), seq(-60, 60, by = 30)) r <- tissot(xy, "+proj=robin") ii <- indicatrix(r) ## Uniform fill plot(ii, scale = 6e5, add = FALSE) tissot_map() ## Colour by areal distortion plot(ii, scale = 6e5, add = FALSE, fill.by = "scale_area") tissot_map()xy <- expand.grid(seq(-150, 150, by = 30), seq(-60, 60, by = 30)) r <- tissot(xy, "+proj=robin") ii <- indicatrix(r) ## Uniform fill plot(ii, scale = 6e5, add = FALSE) tissot_map() ## Colour by areal distortion plot(ii, scale = 6e5, add = FALSE, fill.by = "scale_area") tissot_map()
Generate ellipse coordinates for an indicatrix
ti_ellipse(x, scale = 1e+05, n = 72, ...)ti_ellipse(x, scale = 1e+05, n = 72, ...)
x |
an indicatrix object |
scale |
scaling factor for the ellipse size in projected units |
n |
number of points on the ellipse |
... |
ignored |
A two-column matrix of projected coordinates tracing the ellipse
indicatrix(), plot.indicatrix()
ii <- indicatrix(c(0, 45), "+proj=robin") ell <- ti_ellipse(ii[[1]], scale = 5e5) head(ell) ## draw on a map tissot_map(target = "+proj=robin", add = FALSE) lines(ell)ii <- indicatrix(c(0, 45), "+proj=robin") ell <- ti_ellipse(ii[[1]], scale = 5e5) head(ell) ## draw on a map tissot_map(target = "+proj=robin", add = FALSE) lines(ell)
Compute the Tissot indicatrix at given longitude/latitude locations for a map projection. Returns scale factors, angular deformation, convergence, and related distortion properties.
tissot( x, target, ..., source = "EPSG:4326", method = c("proj", "finitediff"), A = 6378137, f.inv = 298.257223563, dx = 1e-04 )tissot( x, target, ..., source = "EPSG:4326", method = c("proj", "finitediff"), A = 6378137, f.inv = 298.257223563, dx = 1e-04 )
x |
input coordinates — any xy-ish object: a two-column matrix,
data.frame, tibble, list with |
target |
target projection CRS string (required) |
... |
ignored |
source |
source CRS (default |
method |
computation method: |
A |
semi-major axis of the ellipsoid (default WGS84; |
f.inv |
inverse flattening (default WGS84; |
dx |
finite difference step in degrees (default 1e-4; |
By default (method = "proj") distortion measures are computed via
PROJ::proj_factors(), which calls the PROJ C library directly and is
more accurate than finite differences (particularly for tabular or
piecewise-defined projections). The original finite-difference path
based on Snyder (1987) — inspired by Bill Huber's formulation at
https://gis.stackexchange.com/a/5075/482 — is preserved as
method = "finitediff".
x is assumed to contain longitude/latitude values; the default source CRS
is "EPSG:4326". Set source for a different geographic CRS.
A tissot_tbl tibble with columns: x (lon), y (lat),
dx_dlam, dy_dlam, dx_dphi, dy_dphi, scale_h, scale_k,
scale_omega, scale_a, scale_b, scale_area,
angle_deformation, convergence. The source and target CRS
strings are stored as attributes.
indicatrix(), tissot_map(), tissot_unproject(),
PROJ::proj_factors(), PROJ::proj_trans()
tissot(c(0, 45), "+proj=robin") tissot(cbind(seq(-180, 180, by = 30), 0), "+proj=robin") ## compare methods tissot(c(0, 45), "+proj=robin", method = "proj") tissot(c(0, 45), "+proj=robin", method = "finitediff")tissot(c(0, 45), "+proj=robin") tissot(cbind(seq(-180, 180, by = 30), 0), "+proj=robin") ## compare methods tissot(c(0, 45), "+proj=robin", method = "proj") tissot(c(0, 45), "+proj=robin", method = "finitediff")
Deprecated. Prefer passing target explicitly.
These functions access the package-level projection state. Prefer passing
target explicitly to tissot_map() and tissot_abline(), or let
plot.indicatrix_list() or image.tissot_raster() set it automatically.
tissot_get_proj() tissot_set_proj(target)tissot_get_proj() tissot_set_proj(target)
target |
projection CRS string |
tissot_get_proj() returns the current projection string or NULL
tissot_set_proj("+proj=robin") tissot_get_proj() tissot_set_proj(NULL) # resettissot_set_proj("+proj=robin") tissot_get_proj() tissot_set_proj(NULL) # reset
tissot_map() draws the bundled world coastline, projected if a
projection is current. The projection is determined in this order:
An explicit target argument
The projection recorded by the most recent call to
plot.indicatrix_list(), image.tissot_raster(), or
tissot_raster()
tissot_map(..., target = NULL, add = TRUE) tissot_abline(x, y = NULL, ..., source = "EPSG:4326", target = NULL)tissot_map(..., target = NULL, add = TRUE) tissot_abline(x, y = NULL, ..., source = "EPSG:4326", target = NULL)
... |
graphical parameters passed to |
target |
target CRS. If |
add |
logical; add to existing plot (default |
x |
longitude values (or any xy-ish input; see |
y |
latitude values (ignored if |
source |
source CRS for the coordinates (default |
tissot_abline() draws vertical and horizontal reference lines at a
given longitude/latitude in projected coordinates.
tissot_map() invisibly returns the (projected) world coastline matrix
tissot_abline() is called for its side effect
plot.indicatrix_list(), tissot()
r <- tissot(cbind(seq(-150, 150, by = 30), 0), "+proj=robin") ii <- indicatrix(r) plot(ii, scale = 6e5, add = FALSE) tissot_map()r <- tissot(cbind(seq(-150, 150, by = 30), 0), "+proj=robin") ii <- indicatrix(r) plot(ii, scale = 6e5, add = FALSE) tissot_map()
Create a raster-like grid of Tissot distortion metrics in the target
projection. The grid is laid out in projected coordinates, back-projected
to longitude/latitude, and passed through tissot(). Returns a list of
matrices compatible with image() and easily converted to raster objects
(e.g. via terra::rast()).
tissot_raster( target, extent = NULL, radius = 2e+07, nx = 100L, ny = NULL, metrics = c("scale_area", "angle_deformation", "scale_h", "scale_k"), ..., source = "EPSG:4326" ) ## S3 method for class 'tissot_raster' image(x, metric = NULL, col = NULL, asp = 1, main = NULL, ...) ## S3 method for class 'tissot_raster' plot(x, ...)tissot_raster( target, extent = NULL, radius = 2e+07, nx = 100L, ny = NULL, metrics = c("scale_area", "angle_deformation", "scale_h", "scale_k"), ..., source = "EPSG:4326" ) ## S3 method for class 'tissot_raster' image(x, metric = NULL, col = NULL, asp = 1, main = NULL, ...) ## S3 method for class 'tissot_raster' plot(x, ...)
target |
target projection CRS string |
extent |
numeric length 4: |
radius |
maximum half-width of the grid extent in projected units
(default |
nx |
integer; number of grid cells in x (default 100) |
ny |
integer; number of grid cells in y. If |
metrics |
character vector of column names from |
... |
passed to |
source |
source CRS (default |
x |
a |
metric |
which metric to plot (default: first available) |
col |
colour palette (default: |
asp |
aspect ratio (default 1) |
main |
plot title (default: metric name) |
If extent is NULL, the function samples a dense set of boundary
coordinates spanning (longitude/latitude), projects them into the target CRS, and uses
the resulting bounding range (padded by 5%) as the grid extent.
A tissot_raster list with components:
xnumeric vector of x (easting) cell centers
ynumeric vector of y (northing) cell centers
znamed list of ny x nx matrices, one per metric
targetthe target CRS string
extentthe c(xmin, xmax, ymin, ymax) used
The object has an image() method for quick plotting.
tissot(), plot.tissot_raster(), image.tissot_raster()
tr <- tissot_raster("+proj=robin", nx = 60) image(tr) ## specific metric image(tr, metric = "angle_deformation") ## with coastline overlay image(tr) tissot_map() ## polar stereo — use radius to cap extent tp <- tissot_raster("+proj=stere +lat_0=-90", radius = 5e6, nx = 60) image(tp) tissot_map()tr <- tissot_raster("+proj=robin", nx = 60) image(tr) ## specific metric image(tr, metric = "angle_deformation") ## with coastline overlay image(tr) tissot_map() ## polar stereo — use radius to cap extent tp <- tissot_raster("+proj=stere +lat_0=-90", radius = 5e6, nx = 60) image(tp) tissot_map()
A convenience wrapper around PROJ::proj_trans() for converting
projected coordinates to geographic. Useful for generating regular grids
in a projected CRS to feed to tissot(), which requires lon/lat input.
tissot_unproject(x, target = "EPSG:4326", ..., source = NULL)tissot_unproject(x, target = "EPSG:4326", ..., source = NULL)
x |
input coordinates — any xy-ish object: a two-column matrix,
data.frame, tibble, list with |
target |
target CRS string (default |
... |
ignored |
source |
source CRS string (required). Must be a projected CRS. |
A two-column matrix of longitude and latitude values.
## regular grid in UTM zone 55S, unprojected to lon/lat for tissot() xy <- expand.grid(x = seq(4e5, 6e5, length.out = 5), y = seq(5200000, 5400000, length.out = 4)) ll <- tissot_unproject(xy, source = "EPSG:32755") tissot(ll, "+proj=utm +zone=55 +south")## regular grid in UTM zone 55S, unprojected to lon/lat for tissot() xy <- expand.grid(x = seq(4e5, 6e5, length.out = 5), y = seq(5200000, 5400000, length.out = 4)) ll <- tissot_unproject(xy, source = "EPSG:32755") tissot(ll, "+proj=utm +zone=55 +south")
A matrix of longitude/latitude coordinates representing a simplified
world coastline, derived from the maps package. Continent boundaries
are separated by NA rows. Longitudes are constrained to
abs(lon) <= 180.
worldworld
A two-column numeric matrix (longitude, latitude)