## Code

```
library(dplyr)
library(sf)
library(tmap)
library(smoothr) # this is for interpolating along lines
```

figures

code

R

The figure produced below extends to ±89° which makes the point even more clearly, but doesn’t give a nice square Mercator projected map. See also Jason Davies’ page about loxodromes for more on this, including a couple of figures that probably on some level inspired mine. Indeed, Jason Davies’ pages include some gems for exploring the variety of global projections. See, for example, Map Projection Transitions.

This process is pretty complicated to do in *R*, so we need a bunch of helper functions. First, it is the default for `sf`

now, but just to make the point, we intially use S2 mode, so that when we clip data with a hemisphere it does it properly.

Also, define an orthographic projection for the globe view.

Make a hemisphere for the globe view, and apply it to the built-in `World`

dataset that ships with `tmap`

.

Next, convenience ‘helper’ functions for converting coordinate pairs between projections, and converting degrees to radians. By default it will convert longitude-latitude to Mercator.

The next function makes a loxodrome line of equal bearing, starting from 89°S 180°E, and ending when it hits latitude 89°N (latitude range is determined by the `lat`

parameter). The calculations are done in Mercator coordinates, since it is in this projection that a loxodrome is a straight line.

```
get_loxodrome <- function(lat = -89, bearing = 85, n = 100) {
transects <- c() # empty vector for the west to east transects
pt <- c(-180, lat)
# we keep going until we hit the latitude at lat North
while(TRUE && pt[2] < -lat) {
p1 <- x1y1_to_x2y2(pt) # convert to Mercator
# get the point at +180, ie 2 * pi radians in Mercator
p2 <- c(-p1[1],
p1[2] + tan(get_radians(90 - bearing)) * -2 * p1[1])
# the line is a densified version of this
transect <- st_linestring(matrix(c(p1, p2), 2, 2, byrow = TRUE)) %>%
st_sfc(crs = "+proj=merc") %>%
densify(n) # the densification step (provided by smoothr)
transects <- c(transects, transect)
# reset p1 to 'the other side' of the Mercator space i.e. -180
p1 <- p2
pt <- x1y1_to_x2y2(p1, crs1 = "+proj=merc", crs2 = 4326)
pt[1] <- -180
}
# transects need tidying so they extend equally far N and S of equator
# max y coordinate should be the inverse of the minimum y coordinate
ymax <- x1y1_to_x2y2(c(0, -lat))[2]
# apply this limit to the points along the last transect from west to east
n_transects <- length(lines)
# convert the last transect to a set of points, to apply this limit
pts <- transects[[n_transects]] %>%
st_cast("MULTIPOINT") %>%
st_coordinates()
pts <- pts[pts[, 2] <= ymax, 1:2]
# and then convert back to a linestring
transects[[n_transects]] <- pts %>%
matrix(ncol = 2) %>%
st_linestring()
# finally convert to lon-lat i.e. EPSG 4326
transects %>% st_sfc(crs = "+proj=merc") %>%
st_as_sf() %>%
st_transform(4326)
}
```

Make a loxodrome and clip the world to the chosen latitude limits. To apply rectangular projection based limits we have to switch to planar geometry in `sf`

. We have to do this because we can’t show the whole world in Mercator…

```
lox <- get_loxodrome()
lox_o <- lox %>%
st_intersection(hemisphere)
sf_use_s2(FALSE)
mercator_limits <- st_polygon(list(
matrix(c(-180, -89, 180, -89, 180, 89, -180, 89, -180, -89),
ncol = 2, byrow = TRUE))) %>%
st_sfc(crs = 4326) %>%
st_as_sf()
world_m <- World %>%
st_intersection(mercator_limits) %>%
st_transform("+proj=merc") %>%
filter(st_is_empty(geometry) == FALSE)
lox_m <- lox %>%
st_intersection(mercator_limits)
```

And finally make the maps. (In the version in the book I handmade the graticule, not realising that there was a `tm_graticules`

function in `tmap`

…).

```
m1 <- tm_shape(world_o) +
tm_fill(col = "lightgray") +
tm_graticules(x = seq(-180, 165, 15), y = seq(-75, 75, 15),
col = "gray", lwd = 0.5, labels.show = FALSE) +
tm_shape(lox_o) +
tm_lines(col = "black") +
tm_layout(frame = FALSE)
m2 <- tm_shape(world_m, bbox = mercator_limits) +
tm_fill(col = "lightgray") +
tm_graticules(x = seq(-180, 165, 15), y = seq(-75, 75, 15),
col = "gray", lwd = 0.5, labels.show = FALSE) +
tm_shape(lox_m) +
tm_lines(col = "black") +
tm_layout(frame = FALSE)
tmap_arrange(m1, m2, ncol = 2)
```

```
# License (MIT)
#
# Copyright (c) 2023 David O'Sullivan
#
# Permission is hereby granted, free of charge, to any person
# obtaining a copy of this software and associated documentation
# files (the "Software"), to deal in the Software without restriction,
# including without limitation the rights to use, copy, modify, merge,
# publish, distribute, sublicense, and/or sell copies of the Software,
# and to permit persons to whom the Software is furnished to do so,
# subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included
# in all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
# FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
# DEALINGS IN THE SOFTWARE.
```

© 2023 David O’Sullivan