We aim to visualize the creation of map projections in an interactive plot.

First, extract some map data and prepare.

library(maptools)
library(rgl)
library(rglwidget)
library(rgdal)

data(wrld_simpl)

## raw coordinates from maptools
ll <- coordinates(as(as(wrld_simpl, "SpatialLines"), "SpatialPoints"))

## reduce the input map to the south
maxlat <- -10
llsub <- ll[,2] < maxlat

## use the PROJ.4 sphere
a <- 6370997
## function to produce xyz from longitude, latitude, height
## (spherical)
llh2xyz <- function(lonlatheight, rad = 6370997, exag = 1) {
  d2r <- pi / 180.0
  cosLat = cos(lonlatheight[,2] * d2r)
  sinLat = sin(lonlatheight[,2] * d2r)
  cosLon = cos(lonlatheight[,1] * d2r)
  sinLon = sin(lonlatheight[,1] * d2r)
  x = rad * cosLat * sinLon
  y = rad * cosLat * cosLon
  z = (lonlatheight[,3] * exag) + rad * sinLat
  cbind(x, y,-z)
}

xyz <- llh2xyz(cbind(ll, 0), rad = a)

Define three projections, these are Polar Stereographic, Gnomonic and Orthographic all on the South Pole.

Polar Stereographic

family <- c("stere", "gnom", "ortho")
proj <- sprintf("+proj=%s +lon_0=0 +lat_0=-90 +ellps=sphere", family)

jj <- 1

pxy <- project(ll, proj[jj])
  
  ## these are the projected map points on the plane
  pxyz <- cbind(pxy, a)
  
  open3d()
## wgl 
##   1
  ## plot 
  bg3d(bg = "black")
  plot3d(xyz, col = "dodgerblue", axes = FALSE)
  points3d(pxyz[llsub, ], col = "#6AB787FF")
  
  ## rays from the projection point
  ptz2 <- cbind(project(ll[llsub, ][sample(sum(llsub), 100), ], proj[jj]), a)
 
 for (i in seq_len(nrow(ptz2))) {
    origin <- switch(family[jj], 
                    stere = cbind(0, 0, -a), 
                    gnom = cbind(0, 0, 0), 
                    ortho = cbind(ptz2[i, 1], ptz2[i, 2], 0))
    lines3d(rbind(origin, ptz2[i,,drop = FALSE]), color = "grey", lwd =1)
  }
  subid <- currentSubscene3d()
  rglwidget(elementId=family[jj], width = 700, height = 700)

rgl.close()

Gnomonic

jj <- 2

pxy <- project(ll, proj[jj])
## Warning in project(ll, proj[jj]): 19524 projected point(s) not finite
  ## these are the projected map points on the plane
  pxyz <- cbind(pxy, a)
  
  open3d()
## wgl 
##   2
  ## plot 
  bg3d(bg = "black")
  plot3d(xyz, col = "dodgerblue", axes = FALSE)
  points3d(pxyz[llsub, ], col = "#6AB787FF")
  
  ## rays from the projection point
  ptz2 <- cbind(project(ll[llsub, ][sample(sum(llsub), 100), ], proj[jj]), a)
 
 for (i in seq_len(nrow(ptz2))) {
    origin <- switch(family[jj], 
                    stere = cbind(0, 0, -a), 
                    gnom = cbind(0, 0, 0), 
                    ortho = cbind(ptz2[i, 1], ptz2[i, 2], 0))
    lines3d(rbind(origin, ptz2[i,,drop = FALSE]), color = "grey", lwd =1)
  }
  subid <- currentSubscene3d()
  rglwidget(elementId=family[jj], width = 700, height = 700)

rgl.close()

Orthographic

jj <- 3

pxy <- project(ll, proj[jj])
## Warning in project(ll, proj[jj]): 19524 projected point(s) not finite
  ## these are the projected map points on the plane
  pxyz <- cbind(pxy, a)
  
  open3d()
## wgl 
##   3
  ## plot 
  bg3d(bg = "black")
  plot3d(xyz, col = "dodgerblue", axes = FALSE)
  points3d(pxyz[llsub, ], col = "#6AB787FF")
  
  ## rays from the projection point
  ptz2 <- cbind(project(ll[llsub, ][sample(sum(llsub), 100), ], proj[jj]), a)
 
 for (i in seq_len(nrow(ptz2))) {
    origin <- switch(family[jj], 
                    stere = cbind(0, 0, -a), 
                    gnom = cbind(0, 0, 0), 
                    ortho = cbind(ptz2[i, 1], ptz2[i, 2], 0))
    lines3d(rbind(origin, ptz2[i,,drop = FALSE]), color = "grey", lwd =1)
  }
  subid <- currentSubscene3d()
  rglwidget(elementId=family[jj], width = 700, height = 700)

rgl.close()