Figure 3.8 A loxodrome on the sphere and projected

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.

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

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.

sf_use_s2(TRUE)

Also, define an orthographic projection for the globe view.

ortho_proj <- "+proj=ortho lon_0=0 lat_0=40"

Make a hemisphere

Make a hemisphere for the globe view, and apply it to the built-in World dataset that ships with tmap.

hemisphere <- st_point(c(0, 0)) %>%
  st_buffer(6356752) %>%
  st_sfc(crs = ortho_proj) %>%
  densify() %>%
  st_transform(4326)

data("World")
world_o <- World %>%
  st_intersection(hemisphere) %>%
  st_transform(ortho_proj) %>%
  filter(st_is_empty(geometry) == FALSE)

Helper functions for coordinate transformations

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

x1y1_to_x2y2 <- function(coords, crs1 = 4326, crs2 = "+proj=merc") {
  coords %>% st_point() %>%
    st_sfc(crs = crs1) %>%
    st_transform(crs2) %>%
    st_coordinates() %>%
    c()
}

get_radians <- function(d) {
  d * pi / 180
} 

Making a loxodrome

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)
}

Now make the figure!

Assemble the layers

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)

Put them together

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