# Chapter 2 Coordinates

Data are not just numbers, they are numbers with a context”; “In data analysis, context provides meaning(Cobb and Moore 1997)

Before we can try to understand geometries like points, lines, polygons, coverage and grids, it is useful to review coordinate systems so that we have an idea what exactly coordinates of a point reflect. For spatial data, the location of observations are characterized by coordinates, and coordinates are defined in a coordinate system. Different coordinate systems can be used for this, and the most important difference is whether coordinates are defined over a 2-dimensional or 3-dimensional space referenced to orthogonal axes (Cartesian coordinates), or using distance and directions (polar coordinates, spherical and ellipsoidal coordinates). Besides a location of observation, all observations are associated with time of observation, and so time coordinate systems are also briefly discussed. First we will briefly review quantities, to learn what units and datum are.

## 2.1 Quantities, units, datum

The VIM (“International Vocabulary of Metrology”, BIPM et al. (2012)) defines a quantity as a “property of a phenomenon, body, or substance, where the property has a magnitude that can be expressed as a number and a reference”, where “[a] reference can be a measurement unit, a measurement procedure, a reference material, or a combination of such.”

Although one could argue whether all data is constituted of quantities, there is no need to argue that proper data handling requires that numbers (or symbols) are accompanied by information on what they mean, in particular what they refer to.

A measurement system consist of base units for base quantities, and derived units for derived quantities. For instance, the SI system of units (Bureau International des Poids et Mesures 2006) consist of the seven base units length (metre, m), mass (kilogram, kg), time (second, s), electric current (ampere, A), thermodynamic temperature (Kelvin, K), amount of substance (mole, mol), and luminous intensity (candela, cd). Derived units are composed of products of integer powers of base units; examples are speed ($$\mbox{m}~\mbox{s}^{-1}$$), density ($$\mbox{kg}~\mbox{m}^{-3}$$) and area ($$\mbox{m}^2$$).

The special case of unitless measures can refer to either cases where units cancel out (e.g. mass fraction: kg/kg, or angle measured in rad: m/m) or to cases where objects or events were counted (e.g. 5 apples). Adding an angle to a count of apples would not make sense; adding 5 apples to 3 oranges may make sense if the result is reinterpreted as a superclass, e.g. as pieces of fruit. Many data variables have units that are not expressible as SI base units or derived units. Hand (2004) discusses many such measurement scales, e.g. those used to measure intelligence in social sciences, in the context of measurement units.

For many quantities, the natural origin of values is zero. This works for amounts, where differences between amounts result in meaningful negative values. For locations and times, differences have a natural zero interpretation: distance and duration. Absolute location (position) and time need a fixed origin, from which we can meaningfully measure other absolute space-time points: we call this a datum. For space, a datum involves more than one dimension. The combination of a datum and a measurement unit (scale) is a reference system.

We will now elaborate how spatial locations can be expressed as either ellipsoidal or Cartesian coordinates. The next sections will deal with temporal and spatial reference systems, and how they are handled in R.

## 2.2 Ellipsoidal coordinates

par(mar = rep(0,4))
plot(3, 4, xlim = c(-6,6), ylim = c(-6,6), asp = 1)
axis(1, pos = 0, at = 0:6)
axis(2, pos = 0, at = -6:6)
xd = seq(-5, 5, by = .1)
lines(xd, sqrt(25 - xd^2), col = 'grey')
lines(xd, -sqrt(25 - xd^2), col = 'grey')
arrows(0, 0, 3, 4, col = 'red', length = .15, angle = 20)
text(1.5, 2.7, label = "r", col = 'red')
xd = seq(3/5, 1, by = .1)
lines(xd, sqrt(1 - xd^2), col = 'red')
text(1.2, 0.5, label = parse(text = "phi"), col = 'red')
lines(c(3,3), c(0,4), lty = 2, col = 'blue')
lines(c(0,3), c(4,4), lty = 2, col = 'blue')
text(3.3, 0.3, label = "x", col = 'blue')
text(0.3, 4.3, label = "y", col = 'blue')

Figure 2.1 shows both polar and Cartesian coordinates for a two-dimensional situation. In Cartesian coordinates, the point shown is $$(x,y) = (3,4)$$, for polar coordinates it is $$(r,\phi) = (5, \mbox{arctan}(4/3))$$, where $$\mbox{arctan}(4/3)$$ is approximately $$0.93$$ radians, or $$53^{\circ}$$. Note that $$x$$, $$y$$ and $$r$$ all have length units, where $$\phi$$ is an angle (a unitless length/length ratio). Converting back and forth between Cartesian and polar coordinates is trivial, as $x = r~\mbox{cos} \phi,$ $y = r~\mbox{sin} \phi,$ $r = \sqrt{x^2 + y^2}, \ \mbox{and}$ $\phi = \mbox{atan2}(y, x)$ where $$\mbox{atan2}$$ is used in favor of $$\mbox{atan}(y/x)$$ to take care of the right quadrant.

### 2.2.1 Ellipsoidal coordinates

In three dimensions, where Cartesian coordinates are expressed as $$(x,y,z)$$, spherical coordinates are the three-dimensional equivalent of polar coordinates and can be expressed as $$(r,\lambda,\phi)$$, where:

• $$r$$ is the radius of the sphere,
• $$\lambda$$ is the longitude, measured in the $$(x,y)$$ plane counter-clockwise from positive $$x$$, and
• $$\phi$$ is the latitude, the angle between the vector and the $$(x,y)$$ plane.

Figure 2.2 illustrates Cartesian geocentric and ellipsoidal coordinates.

suppressPackageStartupMessages(library(sf))
e = cbind(-90:90,0) # equator
f1 = rbind(cbind(0, -90:90)) # 0/antimerid.
f2 = rbind(cbind(90, -90:90), cbind(270, 90:-90))# +/- 90
eq = st_sfc(st_linestring(e), st_linestring(f1), st_linestring(f2), crs=4326)

geoc = st_transform(eq, "+proj=geocent")
cc = rbind(geoc[[1]], NA, geoc[[2]], NA, geoc[[3]])
from3d = function(x, offset, maxz, minz) {
x = x[,c(2,3,1)] + offset # move to y right, x up, z backw
x[,2] = x[,2] - maxz      # shift y to left
d = maxz
z = x[,3] - minz + offset
x[,1] = x[,1] * (d/z)
x[,2] = x[,2] * (d/z)
x[,1:2]
}
maxz = max(cc[,3], na.rm = TRUE)
minz = min(cc[,3], na.rm = TRUE)
offset = 3e7
circ = from3d(cc, offset, maxz, minz)
mx = max(cc, na.rm = TRUE) * 1.1
x = rbind(c(0, 0, 0), c(mx, 0, 0))
y = rbind(c(0, 0, 0), c(0, mx, 0))
z = rbind(c(0, 0, 0), c(0, 0, mx))
ll = rbind(x, NA, y, NA, z)
l0 =  from3d(ll, offset, maxz, minz)
mx = max(cc, na.rm = TRUE) * 1.2
x = rbind(c(0, 0, 0), c(mx, 0, 0))
y = rbind(c(0, 0, 0), c(0, mx, 0))
z = rbind(c(0, 0, 0), c(0, 0, mx))
ll = rbind(x, NA, y, NA, z)
l =  from3d(ll, offset, maxz, minz)

par(mfrow = c(1, 2))
par(mar = rep(0,4))
plot.new()
plot.window(xlim = c(min(circ[,1],na.rm = TRUE), 3607103*1.02),
ylim = c(min(circ[,2],na.rm = TRUE), 2873898*1.1), asp = 1)
lines(circ)
lines(l0)
text(l[c(2,5,8),], labels = c("x", "y", "z"), col = 'red')
p = st_as_sfc("POINT(60 47)", crs = 4326) %>% st_transform("+proj=geocent")
p = p[[1]]
pts = rbind(c(0,0,0), c(p[1],0,0), c(p[1],p[2],0), c(p[1],p[2],p[2]))
ptsl = from3d(pts, offset, maxz, minz)
lines(ptsl, col = 'blue', lty = 2, lwd = 2)
points(ptsl[4,1], ptsl[4,2], col = 'blue', cex = 1, pch = 16)

plot.new()
plot.window(xlim = c(min(circ[,1],na.rm = TRUE), 3607103*1.02),
ylim = c(min(circ[,2],na.rm = TRUE), 2873898*1.1), asp = 1)
lines(circ)

p = st_as_sfc("POINT(60 47)", crs = 4326) %>% st_transform("+proj=geocent")
p = p[[1]]
pts = rbind(c(0,0,0), c(p[1],p[2],p[3]))
pt =  from3d(pts, offset, maxz, minz)
lines(pt)
points(pt[2,1], pt[2,2], col = 'blue', cex = 1, pch = 16)

p0 = st_as_sfc("POINT(60 0)", crs = 4326) %>% st_transform("+proj=geocent")
p0 = p0[[1]]
pts = rbind(c(0,0,0), c(p0[1],p0[2],p0[3]))
pt =  from3d(pts, offset, maxz, minz)
lines(pt)

p0 = st_as_sfc("POINT(0 0)", crs = 4326) %>% st_transform("+proj=geocent")
p0 = p0[[1]]
pts = rbind(c(0,0,0), c(p0[1],p0[2],p0[3]))
pt =  from3d(pts, offset, maxz, minz)
lines(pt)

p0 = st_as_sfc("POINT(0 90)", crs = 4326) %>% st_transform("+proj=geocent")
p0 = p0[[1]]
pts = rbind(c(0,0,0), c(p0[1],p0[2],p0[3]))
pt =  from3d(pts, offset, maxz, minz)
lines(pt, lty = 2)

p0 = st_as_sfc("POINT(90 0)", crs = 4326) %>% st_transform("+proj=geocent")
p0 = p0[[1]]
pts = rbind(c(0,0,0), c(p0[1],p0[2],p0[3]))
pt =  from3d(pts, offset, maxz, minz)
lines(pt, lty = 2)

f1 = rbind(cbind(0:60, 0))
arc = st_sfc(st_linestring(f1), crs=4326)
geoc = st_transform(arc, "+proj=geocent")
cc = rbind(geoc[[1]])
circ = from3d(cc, offset, maxz, minz)
lines(circ, col = 'red', lwd = 2, lty = 2)

f1 = rbind(cbind(60, 0:47))
arc = st_sfc(st_linestring(f1), crs=4326)
geoc = st_transform(arc, "+proj=geocent")
cc = rbind(geoc[[1]])
circ = from3d(cc, offset, maxz, minz)
lines(circ, col = 'blue', lwd = 2, lty = 2)

text(pt[1,1]+100000, pt[1,2]+50000, labels = expression(phi), col = 'blue') # lat
text(pt[1,1]+20000, pt[1,2]-50000, labels = expression(lambda), col = 'red') # lng

$$\lambda$$ typically varies between $$-180^{\circ}$$ and $$180^{\circ}$$ (or alternatively from $$0^{\circ}$$ to $$360^{\circ}$$), $$\phi$$ from $$-90^{\circ}$$ to $$90^{\circ}$$. When we are only interested in points on a sphere with given radius, we can drop $$r$$: $$(\lambda,\phi)$$ now suffice to identify any point.

It should be noted that this is just a definition, one could for instance also choose to measure polar angle, the angle between the vector and $$z$$, instead of latitude. There is also a long tradition of specifying points as $$(\phi,\lambda)$$ but throughout this book we will stick to longitude-latitude, $$(\lambda,\phi)$$. The point denoted in figure 2.2 has $$(\lambda,\phi)$$ or ellipsoidal coordinates with values

p = st_as_sfc("POINT(60 47)", crs = 4326)
p[[1]]
# POINT (60 47)

with angles measured in degrees, and geocentric coordinates

p = st_as_sfc("POINT(60 47)", crs = 4326) %>% st_transform("+proj=geocent")
p[[1]]
# POINT Z (2178844 3773868 4641765)

with unit metres.

For points on an ellipse, there are two ways in which angle can be expressed (figure 2.3): measured from the center of the ellipse ($$\psi$$), or measured perpendicular to the tangent on the ellipse at the target point ($$\phi$$).

par(mar = rep(0,4))
x = 4
y = 5/8 * sqrt(48)
plot(x, y, xlim = c(-6,6), ylim = c(-8,8), asp = 1)
axis(1, pos = 0, at = 0:9)
axis(2, pos = 0, at = -5:5)
xd = seq(-8, 8, by = .1)
lines(xd, 5/8 * sqrt(64 - xd^2), col = 'grey')
lines(xd, 5/8 * -sqrt(64 - xd^2), col = 'grey')
arrows(0, 0, x, y, col = 'red', length = .15, angle = 20)
b = (x * 25) / (-y * 64)
a = y - x * b
abline(a, b, col = 'grey')
b = -1/b
x0 = x - y / b
arrows(x0, 0, x, y, col = 'blue', length = .15, angle = 20)
text(1.2, 0.5, label = parse(text = "psi"), col = 'red')
text(3, 0.5, label = parse(text = "phi"), col = 'blue')

The most commonly used parametric model for the Earth is an ellipsoid of revolution, an ellipsoid with two equal semi-axes (Iliffe and Lott 2008). In effect, this is a flattened sphere (or spheroid): the distance between the poles is (slightly: about 0.33%) smaller than the distance between two opposite points on the equator. Under this model, longitude is always measured along a circle (as in figure 2.2), and latitude along an ellipse (as in figure 2.3). If we think of figure 2.3 as a cross section of the Earth passing through the poles, the geodetic latitude measure $$\phi$$ is the one used when no further specification is given. The latitude measure $$\psi$$ is called the geocentric latitude.

In addition to longitude and latitude we can add altitude or elevation to define points that are not on the ellipsoid, and obtain a three dimensional space again. When defining altitude, we need to choose:

• where zero altitude is: on the ellipsoid, or relative to the surface approximating mean sea level (the geoid)?
• which direction is positive, and
• which direction is “straight up”: perpendicular to the ellipsoid surface, or in the direction perpendicular to the surface of the geoid?

All these choices may matter, depending on the application area and required measurement accuracies.

The shape of the Earth is not a perfect ellipsoid. As a consequence, several ellipsoids with different shape parameters and bound to the Earth in different ways are being used. Such ellipsoids are called datums, and are briefly discussed in section 2.3, along with coordinate reference systems.

### 2.2.2 Projected coordinates, distances

Because paper maps and computer screens are much more abundant and practical than globes, most of the time we look at spatial data we see it projected: drawn on a flat, two-dimensional surface. Computing the locations in a two-dimensional space means that we work with projected coordinates. Projecting ellipsoidal coordinates means that shapes, directions, areas, or even all three, are distorted (Iliffe and Lott 2008).

Distances between two points $$p_i$$ and $$p_j$$ in Cartesian coordinates are computed as Euclidean distances, in two dimensions by $d_{ij} = \sqrt{(x_i-x_j)^2+(y_i-y_j)^2}$ with $$p_i = (x_i,y_i)$$ and in three dimensions by $d_{ij} = \sqrt{(x_i-x_j)^2+(y_i-y_j)^2+(z_i-z_j)^2}$ with $$p_i = (x_i,y_i,z_i).$$ These distances represent the length of a straight line between two points $$i$$ and $$j$$.

For two points on a circle, the length of the arc of two points $$c_1 = (r,{\phi}_i)$$ and $$c_2 = (r, \phi_2)$$ is $s_{ij}=r~|\phi_1-\phi_2| = r ~\theta$ with $$\theta$$ the angle between $$\phi_1$$ and $$\phi_2$$ in radians. For very small values of $$\theta$$, we will have $$s_{ij} \approx d_{ij}$$, because a small arc segment is nearly straight.

For two points $$p_1 = (\lambda_1,\phi_1)$$ and $$p_2 = (\lambda_2,\phi_2)$$ on a sphere with radius $$r'$$, the great circle distance is the arc length between $$p_1$$ and $$p_2$$ on the circle that passes through $$p_1$$ and $$p_2$$ and has the center of the sphere as its center, and is given by $$s_{12} = r ~ \theta_{12}$$ with $\theta_{12} = \arccos(\sin \phi_1 \cdot \sin \phi_2 + \cos \phi_1 \cdot \cos \phi_2 \cdot \cos(|\lambda_1-\lambda_2|))$ the angle between $$p_1$$ and $$p_2$$, in radians.

Arc distances between two points on a spheroid are more complicated to compute; a good discussion on the topic and an explanation of the method implemented in GeographicLib (part of PROJ) is given in Karney (2013).

To show that these distance measures actually give different values, we computed them for the distance Berlin - Paris. Here, gc_ refers to ellipsoidal and spherical great circle distances, straight_ refers to straight line, Euclidean distances between Cartesian geocentric coordinates associated on the WGS84 ellipse and sphere:

pts = st_sfc(st_point(c(13.4050, 52.5200)), st_point(c(2.3522, 48.8566)), crs = 'EPSG:4326')
s2_orig = sf_use_s2(FALSE)
d1 = c(gc_ellipse = st_distance(pts)[1,2])
sf_use_s2(TRUE)
# or, without using s2, use st_distance(st_transform(pts, "+proj=longlat +ellps=sphere"))
d2 = c(gc_sphere = st_distance(pts)[1,2])
p = st_transform(pts, "+proj=geocent")
d3 = c(straight_ellipse = units::set_units(sqrt(sum(apply(do.call(cbind, p), 1, diff)^2)), m))
p2 = st_transform(pts, "+proj=longlat +ellps=sphere") %>% st_transform("+proj=geocent")
d4 = c(straight_sphere = units::set_units(sqrt(sum(apply(do.call(cbind, p2), 1, diff)^2)), m))
res = c(d1,d3,d2,d4)
# print as km, re-add names:
sf_use_s2(s2_orig) # back to what it was before changing
res %>% units::set_units(km) %>% setNames(names(res)) %>% print(digits = 5)
# Units: [km]
#       gc_ellipse straight_ellipse        gc_sphere  straight_sphere
#           879.70           879.00           877.46           876.77

### 2.2.3 Bounded and unbounded spaces

Two-dimensional and three-dimensional Euclidean spaces ($$R^2$$ and $$R^3$$) are unbounded: every line in this space has infinite length, distances, areas or volumes are unbounded. In contrast, spaces defined on a circle ($$S^1$$) or sphere ($$S^2$$) define a bounded set: there may be infinitely many points but the length and area of the circle and the radius, area and volume of a sphere are bound.

This may sound trivial, but leads to some interesting findings when handling spatial data. A polygon on $$R^2$$ has unambiguously an inside and an outside. On a sphere, $$S^2$$, any polygon divides the sphere in two parts, and which of these two is to be considered inside and which outside is ambiguous and needs to be defined e.g. by the traversal direction. Chapter 4 will further discuss consequences when working with geometries on $$S^2$$.

## 2.3 Coordinate Reference Systems

We follow Lott (2015) when defining the following concepts (italics indicate literal quoting):

• a coordinate system is a set of mathematical rules for specifying how coordinates are to be assigned to points,
• a datum is a parameter or set of parameters that define the position of the origin, the scale, and the orientation of a coordinate system,
• a geodetic datum is a datum describing the relationship of a two- or three-dimensional coordinate system to the Earth, and
• a coordinate reference system is a coordinate system that is related to an object by a datum; for geodetic and vertical datums, the object will be the Earth.

A readable text that further explains these concepts is Iliffe and Lott (2008).

The Earth does not follow a regular shape. The topography of the Earth is of course known to vary strongly, but also the surface formed by constant gravity at mean sea level, the geoid, is irregular. A commonly used model that is fit to the geoid is an ellipsoid of revolution, which is an ellipse with two identical minor axes. Fitting such an ellipsoid to the Earth gives a datum. However, fitting it to different areas, or based on different sets of reference points gives different fits, and hence different datums: a datum can for instance be fixed to a particular tectonic plate (like ETRS89), others can be globally fit (like WGS84). More local fits lead to smaller approximation errors.

The definitions above imply that coordinates in degrees longitude and latitude only have a meaning, i.e. can only be interpreted unambiguously as Earth coordinates, when the datum they are associated with is given.

Note that for projected data, the data that were projected are associated with a reference ellipsoid (datum). Going from one projection to another without changing datum is called coordinate conversion, and passes through the ellipsoidal coordinates associated with the datum involved. This process is lossless and invertible: the parameters and equations associated with a conversion are not empirical. Recomputing coordinates in a new datum is called coordinate transformation, and is approximate: because datums are a result of model fitting, transformations between datums are models too that have been fit; the equations involved are empirical, and multiple transformation paths, based on different model fits and associated with different accuracies, are possible.

Plate tectonics imply that within a global datum, fixed objects may have coordinates that change over time, and that transformations from one datum to another may be time-dependent. Earthquakes are a cause of more local and sudden changes in coordinates.

## 2.4 PROJ and mapping accuracy

Very few living people active in open source geospatial software can remember the time before PROJ. PROJ (Evenden 1990) started in the 1970s as a Fortran project, and was released in 1985 as a C library for cartographic projections. It came with command line tools for direct and inverse projections, and could be linked to software to let it support (re)projection directly. Originally, datums were considered implicit, and no datum transformations were allowed.

In the early 2000s, PROJ was known as PROJ.4, after its never changing major version number. Amongst others motivated by the rise of GPS, the need for datum transformations increased and PROJ.4 was extended with rudimentary datum support. PROJ definitions for coordinate reference systems would look like this:

+proj=utm +zone=33 +datum=WGS84 +units=m +no_defs

where key=value pairs are preceded by a + and separated by a space. This form came to be known as “PROJ.4 string”, since the PROJ project stayed at version 4.x for several decades. Other datums would come with fields like:

+ellps=bessel +towgs84=565.4,50.3,465.6,-0.399,0.344,-1.877,4.072

indicating another ellipse, as well as the seven (or three) parameters for transforming from this ellipse to WGS84 (the “World Geodetic System 1984” global datum once popularized by GPS), effectively defining the datum in terms of a transformation to WGS84.

Along with PROJ.4 came a set of databases with known (registered) projections, from which the best known is the EPSG registry. National mapping agencies would provide (and update over time) their best guesses of +towgs84= parameters for national coordinate reference systems, and distribute it through the EPSG registry, which was part of PROJ distributions. For some transformations, datum grids were available and distributed as part of PROJ.4: such grids are raster maps that provide for every location pre-computed values for the shift in longitude and latitude, or elevation, for a particular datum transformation.

library(stars)
library(rnaturalearth)
countries110 = st_as_sf(countries110)
uk = countries110[countries110\$admin %in% c("United Kingdom"),] %>%
st_geometry()
hook = function() {
plot(uk, border = "orange", col = NA, add = TRUE)
}
plot(r[,,,1:2], axes = TRUE, hook = hook, key.pos = 4)
h = read_stars("data/uk_os_OSGM15_GB.tif")
plot(h, axes = TRUE, reset = FALSE)
plot(uk, border = "orange", col = NA, add = TRUE)

In PROJ.4, every coordinate transformation had to go through a conversion to and from WGS84; even reprojecting data associated with a datum different from WGS84 had to go through a transformation to and from WGS84. The associated errors of up to 100 m were acceptable for mapping purposes for not too small areas, but applications that need high accuracy transformations, e.g. precision agriculture, planning flights of UAV’s, or object tracking are often more demanding in terms of accuracy.

In 2018, after a successful “GDAL Coordinate System Barn Raising” initiative, a number of companies profiting from the open source geospatial software stack supported the development of a more modern, mature coordinate transformation system in PROJ. Over a few years, PROJ.4 evolved through versions 5, 6, 7 and 8 and was hence renamed into PROJ (or PR$$\phi$$J).

The most notable changes include:

• although PROJ.4 strings can still be used to initialize certain coordinate reference systems, they are no longer sufficient to represent all of them; a new format, WKT2 (described in next section) replaces it
• WGS84 as a hub datum is dropped: coordinate transformation no longer need to go through a particular datum
• multiple conversion or transformation paths (so-called pipelines) to go from CRS A to CRS B are possible, and can be reported along with the associated accuracy; PROJ will by default use the most accurate one but user control is possible
• transformation pipelines can chain an arbitrary number of elementary transformation operations, including swapping of axes and unit transformations
• datum grids, of which there are now many more, are no longer distributed with the library but are accessible from a content delivery network (CDN); PROJ allows to enabling and disabling network access to access these grids, and only downloads the section of the grid actually needed, storing it in a cache on the user’s machine for future use
• coordinate transformations receive support for epochs, time-dependent transformations (and hence: four-dimensional coordinates, including the source and target time)
• the set of files with registered coordinate reference systems is handled in an SQLite database
• instead of always handling axis order (longitude, latitude), when the authority defines differently this is now obeyed (with the most notable example: EPSG:4326 defines axis order to be latitude, longitude.)

All these points sound like massive improvements, and accuracies of transformation can be below 1 metre. An interesting point is the last: Where we could safely assume for many decades that spatial data with ellipsoidal coordinates would have axis order (longitude, latitude), this is no longer the case. We will see in section 7.6.6 how to deal with this.

Examples of a horizontal datum grids, downloaded from cdn.proj.org, are shown in figure 2.4 and for a vertical datum grid in figure 2.5. Datum grids may carry per-pixel accuracy values.

## 2.5 WKT-2

Lott (2015) describes a standard for encoding coordinate reference systems, as well as transformations between them using well known text; the standard (and format) is referred to informally as WKT-2. As mentioned above, GDAL and PROJ fully support this encoding. An example of WKT2 for CRS OGC:CRS84 is:

GEOGCRS["WGS 84",
DATUM["World Geodetic System 1984",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]],
ID["EPSG",6326]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8901]],
CS[ellipsoidal,2],
AXIS["longitude",east,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433,
ID["EPSG",9122]]],
AXIS["latitude",north,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433,
ID["EPSG",9122]]]]

This shows a WGS84 ellipsoid, and a coordinate system with the axis order (longitude, latitude) that can be used to replace EPSG:4326 when one wants unambiguously “traditional” (GIS) axis order.

A longer introduction on the history and recent changes in PROJ is given in R. Bivand (2020b), building upon the work of Knudsen and Evers (2017) and Evers and Knudsen (2017).

## 2.6 Exercises

Try to solve the following exercises with R (without loading packages); try to use functions where appropriate:

1. list three geographic measures that do not have a natural zero origin
2. convert the $$(x,y)$$ points $$(10,2)$$, $$(-10,-2)$$, $$(10,-2)$$ and $$(0,10)$$ to polar coordinates
3. convert the polar $$(r,\phi)$$ points $$(10,45^{\circ})$$, $$(0,100^{\circ})$$ and $$(5,359^{\circ})$$ to Cartesian coordinates
4. assuming the Earth is a sphere with a radius of 6371 km, compute for $$(\lambda,\phi)$$ points the great circle distance between $$(10,10)$$ and $$(11,10)$$, between $$(10,80)$$ and $$(11,80)$$, between $$(10,10)$$ and $$(10,11)$$ and between $$(10,80)$$ and $$(10,81)$$ (units: degree). What are the distance units?

### References

BIPM, IEC, ILAC IFCC, IUPAP IUPAC, and OIML ISO. 2012. “The International Vocabulary of Metrology–Basic and General Concepts and Associated Terms (Vim), 3rd Edn. JCGM 200: 2012.” JCGM (Joint Committee for Guides in Metrology). https://www.bipm.org/en/publications/guides/.

Bivand, Roger. 2020b. Why Have Crs, Projections and Transformations Changed? https://rgdal.r-forge.r-project.org/articles/CRS_projections_transformations.html.

Bureau International des Poids et Mesures. 2006. The International System of Units (SI), 8th Edition. Organisation Intergouvernementale de la Convention du Mètre. https://www.bipm.org/en/publications/si-brochure/download.html.

Cobb, George W., and David S. Moore. 1997. “Mathematics, Statistics and Teaching.” The American Mathematical Monthly 104 (9): 801–23. http://www.jstor.org/stable/2975286.

Evenden, Gerald I. 1990. Cartographic Projection Procedures for the UNIX Environment — a User’s Manual. http://download.osgeo.org/proj/OF90-284.pdf.

Evers, Kristian, and Thomas Knudsen. 2017. Transformation Pipelines for Proj.4. https://www.fig.net/resources/proceedings/fig_proceedings/fig2017/papers/iss6b/ISS6B_evers_knudsen_9156.pdf.

Hand, David J. 2004. Measurement: Theory and Practice. A Hodder Arnold Publication.

Iliffe, Jonathan, and Roger Lott. 2008. Datums and Map Projections for Remote Sensing, Gis, and Surveying. Whittles Pub. CRC Press, Scotland, UK.

Karney, Charles FF. 2013. “Algorithms for Geodesics.” Journal of Geodesy 87 (1): 43–55. https://link.springer.com/content/pdf/10.1007/s00190-012-0578-z.pdf.

Knudsen, Thomas, and Kristian Evers. 2017. Transformation Pipelines for Proj.4. https://meetingorganizer.copernicus.org/EGU2017/EGU2017-8050.pdf.

Lott, Roger. 2015. “Geographic Information-Well-Known Text Representation of Coordinate Reference Systems.” Open Geospatial Consortium. http://docs.opengeospatial.org/is/12-063r5/12-063r5.html.