Vous produisez une carte quinconciale Peirce? [fermé]

Pour autant que je sache, ni les outils PROJ4 ni ESRI ne peuvent appliquer la projection quinconciale de Peirce.

Quelqu'un sait-il quelles bibliothèques / logiciels peuvent le gérer?

J'ai récemment remarqué que la lib D3js peut reprojeter à la projection Peirce Quincuncial (avec le plugin geo):

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

Dans R, on peut utiliser cette fonction (copiée ci-dessous) pour transformer chaque coordonnée dans un fichier de formes, puis tracer la carte.

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

Maintenant, comment l'utiliser.

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 est un logiciel commercial qui semble prendre en charge la projection en quinconce de Peirce. (Je ne l'ai pas utilisé moi-même, donc je ne peux pas vérifier comment cela fonctionne.)

Je vois que cette projection est également utilisée pour créer un certain type de photo panoramique . Si vous avez seulement besoin de projeter une image (par opposition aux ensembles de données vectorielles), vous pourrez peut-être trouver une solution de traitement d'image. Par exemple, voici un tutoriel sur la création de panoramas en quinconce Peirce avec des plugins Photoshop, et voici une discussion (avec des scripts) pour appliquer la projection aux images avec MathMap .

L'article Warping Peirce Quincuncial Panoramas de Chamberlain Fong et Brian K. Vogel comprend une implémentation MatLab de leur approche. Il est également orienté image, mais MatLab peut lire des fichiers de formes , donc peut-être qu'une projection vectorielle pourrait être bricolée…