How to make Peirce projection algorithm allow for 360 degrees rotation?

Question

I'm using this algorithm for Peirce world map projection in R. I'm able to do some fine maps, for instance using 28 as the value for the lambda_0 parameter in function toPeirceQuincuncial, since this angle creates less land distortion and breaks no important islands (besides Antarctica, obviously). The algorithm is used like this:

# constants
pi<-acos(-1.0)
twopi<-2.0*pi
halfpi<-0.5*pi
degree<-pi / 180
halfSqrt2<-sqrt(2) / 2
quarterpi<-0.25 * pi
mquarterpi<--0.25 * pi
threequarterpi<-0.75 * pi
mthreequarterpi<--0.75 * pi
radian<-180/pi
sqrt2<-sqrt(2)
sqrt8<-2. * sqrt2
halfSqrt3<-sqrt(3) / 2
PeirceQuincuncialScale<-3.7081493546027438 ;# 2*K(1/2)
PeirceQuincuncialLimit<-1.8540746773013719 ;# K(1/2)

ellFaux<-function(cos_phi,sin_phi,k){
  x<-cos_phi * cos_phi
  y<-1.0 - k * k * sin_phi * sin_phi
  z<-1.0
  rf<-ellRF(x,y,z)
  return(sin_phi * rf)
}

ellRF<-function(x,y,z){
  if (x < 0.0 || y < 0.0 || z < 0.0) {
    print("Negative argument to Carlson's ellRF")
    print("ellRF negArgument")
  }
  delx<-1.0; 
  dely<-1.0; 
  delz<-1.0
  while(abs(delx) > 0.0025 || abs(dely) > 0.0025 || abs(delz) > 0.0025) {
    sx<-sqrt(x)
    sy<-sqrt(y)
    sz<-sqrt(z)
    len<-sx * (sy + sz) + sy * sz
    x<-0.25 * (x + len)
    y<-0.25 * (y + len)
    z<-0.25 * (z + len)
    mean<-(x + y + z) / 3.0
    delx<-(mean - x) / mean
    dely<-(mean - y) / mean
    delz<-(mean - z) / mean
  }
  e2<-delx * dely - delz * delz
  e3<-delx * dely * delz
  return((1.0 + (e2 / 24.0 - 0.1 - 3.0 * e3 / 44.0) * e2+ e3 / 14) / sqrt(mean))
}

toPeirceQuincuncial<-function(lambda,phi,lambda_0=20.0){
  # Convert latitude and longitude to radians relative to the
  # central meridian

  lambda<-lambda - lambda_0 + 180
  if (lambda < 0.0 || lambda > 360.0) {
    lambda<-lambda - 360 * floor(lambda / 360)
  }
  lambda<-(lambda - 180) * degree
  phi<-phi * degree

  # Compute the auxiliary quantities 'm' and 'n'. Set 'm' to match
  # the sign of 'lambda' and 'n' to be positive if |lambda| > pi/2

  cos_phiosqrt2<-halfSqrt2 * cos(phi)
  cos_lambda<-cos(lambda)
  sin_lambda<-sin(lambda)
  cos_a<-cos_phiosqrt2 * (sin_lambda + cos_lambda)
  cos_b<-cos_phiosqrt2 * (sin_lambda - cos_lambda)
  sin_a<-sqrt(1.0 - cos_a * cos_a)
  sin_b<-sqrt(1.0 - cos_b * cos_b)
  cos_a_cos_b<-cos_a * cos_b
  sin_a_sin_b<-sin_a * sin_b
  sin2_m<-1.0 + cos_a_cos_b - sin_a_sin_b
  sin2_n<-1.0 - cos_a_cos_b - sin_a_sin_b
  if (sin2_m < 0.0) {
    sin2_m<-0.0
  }
  sin_m<-sqrt(sin2_m)
  if (sin2_m > 1.0) {
    sin2_m<-1.0
  }
  cos_m<-sqrt(1.0 - sin2_m)
  if (sin_lambda < 0.0) {
    sin_m<--sin_m
  }
  if (sin2_n < 0.0) {
    sin2_n<-0.0
  }
  sin_n<-sqrt(sin2_n)
  if (sin2_n > 1.0) {
    sin2_n<-1.0 
  }
  cos_n<-sqrt(1.0 - sin2_n)
  if (cos_lambda > 0.0) {
    sin_n<--sin_n
  }

  # Compute elliptic integrals to map the disc to the square

  x<-ellFaux(cos_m,sin_m,halfSqrt2)
  y<-ellFaux(cos_n,sin_n,halfSqrt2)

  # Reflect the Southern Hemisphere outward

  if(phi < 0) {
    if (lambda < mthreequarterpi) {
      y<-PeirceQuincuncialScale - y
    } else if (lambda < mquarterpi) {
      x<--PeirceQuincuncialScale - x
    } else if (lambda < quarterpi) {
      y<--PeirceQuincuncialScale - y
    } else if (lambda < threequarterpi) {
      x<-PeirceQuincuncialScale - x
    } else {
      y<-PeirceQuincuncialScale - y
    }
  }

  # Rotate the square by 45 degrees to fit the screen better

  X<-(x - y) * halfSqrt2
  Y<-(x + y) * halfSqrt2
  res<-list(X,Y)
  return(res)
}

library(rgdal)
ang <- 28
p <- readOGR('~/R/shp','TM_WORLD_BORDERS-0.3') # read world shapefile
for (p1 in 1:length(p@polygons)) {
  for (p2 in 1:length(p@polygons[[p1]]@Polygons)) {
    for (p3 in 1:nrow(p@polygons[[p1]]@Polygons[[p2]]@coords)) {
      pos <- toPeirceQuincuncial(p@polygons[[p1]]@Polygons[[p2]]@coords[p3,1],
                                 p@polygons[[p1]]@Polygons[[p2]]@coords[p3,2],ang)
      p@polygons[[p1]]@Polygons[[p2]]@coords[p3,1] <- pos[[1]][1]
      p@polygons[[p1]]@Polygons[[p2]]@coords[p3,2] <- pos[[2]][1]
    }
  }
}
z <- toPeirceQuincuncial(0,-90,ang)[[1]][1]
p@bbox[,1] <- -z
p@bbox[,2] <- z
# plotting the map
par(mar=c(0,0,0,0),bg='#a7cdf2',xaxs='i',yaxs='i')
plot(p,col='gray',lwd=.5)
for (lon in 15*1:24) { # meridians
  pos <- 0
  posAnt <- 0
  for (lat in -90:90) {
    if (length(pos) == 2) {
      posAnt <- pos
    }
    pos <- toPeirceQuincuncial(lon,lat,ang)
    if (length(posAnt) == 2) {
      segments(pos[[1]][1],pos[[2]][1],posAnt[[1]][1],posAnt[[2]][1],col='white',lwd=.5)
    }
  }
}
lats <- 15*1:5
posS <- matrix(0,length(lats),2)
posST <- 0
pos0 <- 0
posN <- matrix(0,length(lats),2)
posNT <- 0
for (lon in 0:360) {
  posAntS <- posS
  posAntST <- posST
  posAnt0 <- pos0
  posAntN <- posN
  posAntNT <- posNT
  pos0 <- unlist(toPeirceQuincuncial(lon,0,ang))
  posST <- unlist(toPeirceQuincuncial(lon,-23.4368,ang))
  posNT <- unlist(toPeirceQuincuncial(lon,23.4368,ang))
  for (i in 1:length(lats)) {
    posS[i,] <- unlist(toPeirceQuincuncial(lon,-lats[i],ang))
    posN[i,] <- unlist(toPeirceQuincuncial(lon,lats[i],ang))
  }
  if (lon > 0) {
    segments(pos0[1],pos0[2],posAnt0[1],posAnt0[2],col='red',lwd=1)
    segments(posNT[1],posNT[2],posAntNT[1],posAntNT[2],col='yellow')
    for (i in 1:length(lats)) {
      segments(posN[i,1],posN[i,2],posAntN[i,1],posAntN[i,2],col='white',lwd=.5)
    }
    if (!(lon %in% round(90*(0:3+.5)+ang))) {
      for (i in 1:length(lats)) {
        segments(posS[i,1],posS[i,2],posAntS[i,1],posAntS[i,2],col='white',lwd=.5)
      }
      segments(posST[1],posST[2],posAntST[1],posAntST[2],col='yellow')
    } else {
      for (i in 1:length(lats)) {
        posS[i,] <- unlist(toPeirceQuincuncial(lon-0.001,-lats[i],ang))
        segments(posS[i,1],posS[i,2],posAntS[i,1],posAntS[i,2],col='white',lwd=.5)
        posS[i,] <- unlist(toPeirceQuincuncial(lon,-lats[i],ang))
      }
      posST <- unlist(toPeirceQuincuncial(lon-0.001,-23.4368,ang))
      segments(posST[1],posST[2],posAntST[1],posAntST[2],col='yellow')
      posST <- unlist(toPeirceQuincuncial(lon,-23.4368,ang))
    }
  }
}

Playing with different values for lambda_0, I've found out that I apparently cannot choose any value I want. It seems that the function will only work with half the possibilities I thought it did.

Numbers indicate values of lambda_0. As you can see, North America moves clockwise from right to left between 80° and 200°, and starts the same movement again between 240° and 40°.

How can I change the algorithm to allow for any angle I want (for instance, North America pointing up)?

Answer 1

I made a simple hack that "solved" the problem. First I've found out which angles were a repetition of the others: between 46 and 225 degrees. Then, for these angles, I just had to flip the bounding box of the shapefile:

z <- toPeirceQuincuncial(0,-90,ang)[[1]][1]
if (ang > 45 & ang < 226) {
  p@bbox[,1] <- z
  p@bbox[,2] <- -z
} else {
  p@bbox[,1] <- -z
  p@bbox[,2] <- z
}

Other thing: I improved the R code of toPeirceQuincuncial, since it was converting coordinate by coordinate. Since R is a vector language, I adapted it to convert a group of coordinates at once, which made the code extremely fast. Now lambda and phi may both be vectors (same size, of course), while lambda_0 is still a single number.

toPeirceQuincuncial2<-function(lambda,phi,lambda_0=20.0){
  # Convert latitude and longitude to radians relative to the
  # central meridian

  lambda<-lambda - lambda_0 + 180
  lambda[which(lambda<0.0 | lambda>360.0)] <-
    lambda[which(lambda<0.0 | lambda>360.0)] - 360*floor(lambda[which(lambda<0.0 | lambda>360.0)]/360)
  lambda<-(lambda - 180) * degree
  phi<-phi * degree

  # Compute the auxiliary quantities 'm' and 'n'. Set 'm' to match
  # the sign of 'lambda' and 'n' to be positive if |lambda| > pi/2

  cos_phiosqrt2<-halfSqrt2 * cos(phi)
  cos_lambda<-cos(lambda)
  sin_lambda<-sin(lambda)
  cos_a<-cos_phiosqrt2 * (sin_lambda + cos_lambda)
  cos_b<-cos_phiosqrt2 * (sin_lambda - cos_lambda)
  sin_a<-sqrt(1.0 - cos_a * cos_a)
  sin_b<-sqrt(1.0 - cos_b * cos_b)
  cos_a_cos_b<-cos_a * cos_b
  sin_a_sin_b<-sin_a * sin_b
  sin2_m<-1.0 + cos_a_cos_b - sin_a_sin_b
  sin2_n<-1.0 - cos_a_cos_b - sin_a_sin_b
  sin2_m[which(sin2_m < 0.0)] <- 0.0
  sin_m<-sqrt(sin2_m)
  sin2_m[which(sin2_m > 1.0)] <- 1.0
  cos_m<-sqrt(1.0 - sin2_m)
  sin_m[which(sin_lambda < 0.0)] <- -sin_m[which(sin_lambda < 0.0)]
  sin2_n[which(sin2_n < 0.0)] <- 0.0
  sin_n<-sqrt(sin2_n)
  sin2_n[which(sin2_n > 1.0)] <- 1.0
  cos_n<-sqrt(1.0 - sin2_n)
  sin_n[which(cos_lambda > 0.0)] <- -sin_n[which(cos_lambda > 0.0)]

  # Compute elliptic integrals to map the disc to the square

  x<-ellFaux(cos_m,sin_m,halfSqrt2)
  y<-ellFaux(cos_n,sin_n,halfSqrt2)

  # Reflect the Southern Hemisphere outward

  y[which(phi<0 & lambda<mthreequarterpi)] <- PeirceQuincuncialScale - y[which(phi<0 & lambda<mthreequarterpi)]
  x[which(phi<0 & lambda>=mthreequarterpi & lambda<mquarterpi)] <- -PeirceQuincuncialScale - x[which(phi<0 & lambda>=mthreequarterpi & lambda<mquarterpi)]
  y[which(phi<0 & lambda>=mquarterpi & lambda<quarterpi)] <- -PeirceQuincuncialScale - y[which(phi<0 & lambda>=mquarterpi & lambda<quarterpi)]
  x[which(phi<0 & lambda>=quarterpi & lambda<threequarterpi)] <- PeirceQuincuncialScale - x[which(phi<0 & lambda>=quarterpi & lambda<threequarterpi)]
  y[which(phi<0 & lambda>=threequarterpi)] <- PeirceQuincuncialScale - y[which(phi<0 & lambda>=threequarterpi)]

  # Rotate the square by 45 degrees to fit the screen better

  X<-(x - y) * halfSqrt2
  Y<-(x + y) * halfSqrt2
  res<-list(X,Y)
  return(res)
}