(02) 坐标系（空间数据科学——R语言应用）

本文由SCY翻译自《Spatial Data Science With Applications in R》第二章

坐标系

“数据不仅仅是数字，它们是带有背景的数字”；“在数据分析中，背景赋予了数据以意义。”[@cobbmoore]

在我们尝试理解点、线、多边形、覆盖层和网格等几何形状之前，回顾坐标系统是有益的，这样我们才能理解一个点的坐标到底反映了什么。对于空间数据而言，观测位置由坐标来表征，而坐标是在坐标系统中定义的。可以使用不同的坐标系统来实现这一点，最重要的区别在于坐标是在二维还是三维空间中定义，这个空间参照正交轴（笛卡尔坐标），或使用距离和方向（极坐标、球面坐标和椭球坐标）。除了观测位置，所有观测都与观测时间相关联，因此时间坐标系统也会简要讨论。首先我们将简要回顾量度，以了解单位和基准是什么。

量度、单位、基准

VIM（“国际计量词汇”，@vim）定义量度为“一个现象、物体或物质的属性，这种属性具有可以用一个数字和参照物表示的大小”，其中“[a]参照物可以是测量单位、测量程序、参照物质或这些的组合。”尽管可以争论所有数据是否都由量度构成，但没有必要争论正确处理数据要求数字（或符号）必须伴随关于它们的含义的信息，尤其是它们指向什么。

一个测量系统由基本单位组成，用于基本量度，以及由衍生单位组成，用于衍生量度。例如，国际单位制[@SI]由七个基本单位组成：长度（米，m）、质量（千克，kg）、时间（秒，s）、电流（安培，A）、热力学温度（开尔文，K）、物质的量（摩尔，mol）和发光强度（坎德拉，cd）。衍生单位由基本单位的整数次幂的乘积组成；例如速度（$\mbox{m}~\mbox{s}^{-1}$）、密度（$\mbox{kg}~\mbox{m}^{-3}$）和面积（$\mbox{m}^2$）。

---

此翻译保留了原文的专业术语和格式，确保其科学性和准确性。如需进一步翻译或有其他问题，请随时告知。

The special case of unitless measures can refer to either cases
where units cancel out (for instance mass fraction: kg/kg, or angle measured
in rad: m/m) or to cases where objects or events were counted
(such as “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 in terms of a superclass, in this case as pieces of fruit.
Many data variables have units that are not expressible as SI base
units or derived units. @hand discusses many such measurement scales,
including those used to measure variables like 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.

\index{datum}

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.

Ellipsoidal coordinates

fig-polar, echo

#| out.width: 60%
#| fig.cap: "Two-dimensional polar (red) and Cartesian (blue) coordinates"
#| code-fold: true
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')

@fig-polar 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, \ \mbox{and}$$
$$r = \sqrt{x^2 + y^2}, \ \ \ \phi = \mbox{atan2}(y, x)$$
where $\mbox{atan2}$ is used in favour of $\mbox{atan}(y/x)$ to take care
of the right quadrant.

Spherical or ellipsoidal coordinates

\index{ellipsoidal coordinates}
\index{coordinates!ellipsoidal}
\index{coordinates!geocentric}
\index{coordinates!Cartesian}

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.

@fig-sphere illustrates Cartesian geocentric and
ellipsoidal coordinates.

fig-sphere, echo

#| fig.cap: "Cartesian geocentric coordinates (left) measure three distances, ellipsoidal coordinates (right) measure two angles, and possibly an ellipsoidal height"
#| code-fold: true
library(sf) |> suppressPackageStartupMessages()
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='OGC:CRS84')

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')
# add POINT(60 47)
p <- st_as_sfc("POINT(60 47)", crs = 'OGC:CRS84') |> 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 = 'OGC:CRS84') |> 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 = 'OGC:CRS84') |> 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 = 'OGC:CRS84') |> 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 = 'OGC:CRS84') |> 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 = 'OGC:CRS84') |> 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='OGC:CRS84')
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='OGC:CRS84')
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.

\newpage
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 @fig-sphere has $(\lambda,\phi)$ or ellipsoidal
coordinates with angular values

echo

#| code-fold: true
#| collapse: false
p <- st_as_sfc("POINT(60 47)", crs = 'OGC:CRS84')
p[[1]]

measured in degrees, and geocentric coordinates

echo

#| code-fold: true
#| collapse: false
p <- st_as_sfc("POINT(60 47)", crs = 'OGC:CRS84') |> st_transform("+proj=geocent")
p[[1]]

measured in metres.

\index{coordinates!units}

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

fig-ellipse, echo

#| out.width: 60%
#| fig.cap: "Angles on an ellipse: geodetic (blue) and geocentric (red) latitude"
#| code-fold: true
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
[@iliffelott]. 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
@fig-sphere), and latitude along an ellipse (as in
@fig-ellipse). If we
think of @fig-ellipse 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.

\index{latitude!geodetic}
\index{latitude!geocentric}
\index{longitude}
\index{ellipsoid of revolution}
\index{altitude}
\index{altitude!direction}
\index{coordinates!altitude}

we can add altitude or elevation to longitude and latitude to
define points that are above or below 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 of gravity, 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 @sec-crs, along
with coordinate reference systems.

Projected coordinates, distances