¿Produciendo el mapa quincuncial de Peirce? [cerrado]

Hasta donde sé, ni las herramientas PROJ4 ni ESRI pueden aplicar la proyección quincuncial de Peirce.

¿Alguien sabe qué bibliotecas / softwares pueden gestionarlo?

Recientemente me di cuenta de que la biblioteca D3js puede reproyectarse a la proyección Quiruncial de Peirce (con el plugin geo):

http://bl.ocks.org/mbostock/4310087

En R, uno puede usar esta función (copiada a continuación) para transformar cada coordenada en un archivo de forma y luego trazar el mapa.

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

Ahora cómo usarlo.

library(rgdal)
p <- readOGR('../shp/ne_110m_admin_0_map_units','ne_110m_admin_0_map_units') # downloaded from https://www.naturalearthdata.com/http//www.naturalearthdata.com/download/110m/cultural/ne_110m_admin_0_map_units.zip
ang <- 28 # the lambda_0 from the Peirce function
# change all coordinates
for (p1 in 1:length(p@polygons)) {
  print(paste0(p1,'/',length(p@polygons)))
  flush.console()
  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]
    }
  }
}
# change the bbox of the SpatialPolygonsDataFrame object (p).
z <- toPeirceQuincuncial(0,-90,ang)[[1]][1]
p@bbox[1,1] <- -z
p@bbox[1,2] <- z
p@bbox[2,1] <- -z
p@bbox[2,2] <- z
# start plotting
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 # parallels
posS <- matrix(0,length(lats),2) # southern parallels
posST <- 0 # southern tropic (Tropic of Capricorn)
pos0 <- 0 # Equator
posN <- matrix(0,length(lats),2) # northern parallels
posNT <- 0 # northern tropic (Tropic of Cancer)
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))
    }
  }
}
dev.print(width=1000,height=1000,'Peirce.png',dev=png)

Mapthematics Geocart es un software comercial que parece apoyar la proyección quincuncial de Peirce. (No lo he usado yo mismo, así que no puedo verificar cómo funciona).

Veo que esta proyección también se usa para crear un cierto tipo de foto panorámica . Si solo necesita proyectar una imagen (a diferencia de los conjuntos de datos vectoriales), puede encontrar una solución de procesamiento de imágenes. Por ejemplo, aquí hay un tutorial sobre cómo crear panoramas quincunciales de Peirce con complementos de Photoshop, y aquí hay una discusión (con scripts) para aplicar la proyección a imágenes con MathMap .

El documento Warping Peirce Quincuncial Panoramas de Chamberlain Fong y Brian K. Vogel incluye una implementación de MatLab de su enfoque. También está orientado a la imagen, pero MatLab puede leer archivos de forma , por lo que tal vez una proyección vectorial podría combinarse ...