search

mdsumner / pproj

Illustrate projections
Other
4 stars 1 forks source link

laea (nearly) #2

Open mdsumner opened 5 years ago

mdsumner commented 5 years ago 
shell <- reproj::reproj(geosphere::randomCoordinates(10000), 
                        "+proj=geocent +datum=WGS84", source = 4326)

zoffs <- 6378137
pts <- geosphere::randomCoordinates(3)
m <- maps::map(plot = FALSE)
pts <- cbind(m$x, m$y)
pts <- pts[!is.na(pts[,1]), ]
pts <- pts[sample(1:nrow(pts), 100), ]

pts <- cbind(sample(1:180, 5), 0)

pts_xyz <- reproj::reproj(pts, "+proj=geocent +datum=WGS84", source = 4326)

laea <- "+proj=laea  +lat_0=-90 +datum=WGS84"
pts_xy <- rgdal::project(pts, laea)

library(rgl)

clear3d()
points3d(shell, col = "firebrick", size = 1)
points3d(cbind(pts_xy, zoffs), col = "green")

for (i in seq_len(nrow(pts_xyz))) {
  line1 <- cbind(pts[i, 1], seq(90,-90, length.out = 36))
  line2 <- cbind(pts[i, 1]-180, seq(90,-90, length.out = 36))
  line <- rbind(line1, line2)

  line_xyz <- reproj::reproj(line, sprintf("+proj=geocent +a=%f", sqrt(sum((c(0, 0, zoffs) - pts_xyz[i, ])^2))), 
                             source = 4326)
  line_xyz[,3] <- line_xyz[,3] + zoffs

  lines3d(line_xyz, alpha = 0.1)

  lines3d(rbind(c(0, 0, 0), pts_xyz[i, ]))
  points3d(pts_xyz[i, , drop = FALSE], size = 8)

}
library(raster)
rr <- raster::raster(raster::extent(-6378137, 6378137, -6378137, 6378137)*4)
nrow(rr) <- 1
ncol(rr) <- 1

plane <- quadmesh::quadmesh(rr)
plane$vb[3, ] <- zoffs
plane$material$color <- rgb(0.2, 0.2, 0.2)
shade3d(plane, alpha = 0.2)
aspect3d("iso"); rglwidget()
mdsumner commented 5 years ago

This:

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)
}
library(sp)
library(rgl)
data("wrld_simpl", package = "maptools")

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

a <- 6378137
xyz <- llh2xyz(cbind(ll, 0), rad = a)
laea <- "+proj=laea +lat_0=-90 +datum=WGS84"
pxy <- reproj::reproj(ll, laea, source = 4326)[,1:2]
## these are the projected map points on the plane
pxyz <- cbind(pxy, a)

clear3d()
## plot 
bg3d(bg = "white")
plot3d(xyz, col = "dodgerblue", axes = FALSE)
points3d(pxyz, col = "#6AB787FF")
aspect3d("iso")
for (i in sample(seq_len(nrow(xyz)), 50)) {
  line <- cbind(ll[i, 1], seq(90,0, length.out = 36))
  dxyz <- c(0, 0, a) - xyz[i, ]
  line_xyz <- llh2xyz(cbind(line, 0), rad = sqrt(sum(dxyz^2)))
  line_xyz[,3] <- line_xyz[,3] + a

  lines3d(line_xyz, alpha = 1, col = "white")

  lines3d(rbind(c(0, 0, a), xyz[i, ]), col = "yellow")
  points3d(xyz[i, , drop = FALSE], size = 8, col = "firebrick")

}

(using PROJ for geocent gives a kind of confusing orientation, for another day)

image