Memproduksi peta quincuncial Peirce? [Tutup]

Sejauh yang saya tahu, baik alat PROJ4 maupun ESRI tidak dapat menerapkan proyeksi quincuncial Peirce.

Adakah yang tahu perpustakaan / perangkat lunak apa yang dapat mengelolanya?

Saya baru-baru ini memperhatikan bahwa lib D3js dapat memproyeksikan kembali ke proyeksi Quincuncial Peirce (dengan plugin geo):

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

Di R, seseorang dapat menggunakan fungsi ini (disalin di bawah) untuk mengubah setiap koordinat dalam sebuah shapefile, dan kemudian memplot peta.

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

Sekarang bagaimana menggunakannya.

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 adalah perangkat lunak komersial yang muncul untuk mendukung proyeksi quincuncial Peirce. (Saya belum menggunakannya sendiri, jadi saya tidak bisa memverifikasi cara kerjanya.)

Saya melihat proyeksi ini juga digunakan untuk membuat foto panorama jenis tertentu . Jika Anda hanya perlu memproyeksikan gambar (tidak seperti dataset vektor), Anda mungkin dapat menemukan solusi pemrosesan gambar. Sebagai contoh, berikut ini adalah tutorial tentang cara membuat panorama quincuncial Peirce dengan plugin Photoshop, dan ini adalah diskusi (dengan skrip) untuk menerapkan proyeksi gambar dengan MathMap .

Makalah Warping Peirce Quincuncial Panoramas oleh Chamberlain Fong dan Brian K. Vogel termasuk implementasi MatLab dari pendekatan mereka. Ini berorientasi pada gambar juga, tetapi MatLab dapat membaca shapefile , jadi mungkin proyeksi vektor dapat dirangkai menjadi satu ...