meshblah.R

globe <- function(x, gproj = "+proj=geocent +ellps=WGS84", p4 = "+init=epsg:4326", ...) {
  haveZ <- "z_" %in% names(x$v)
  
  ## need to handle if we already have a "z_"
  if (haveZ) {
    ll <- as.matrix(x[, c("x_", "y_", "z_")])
    
  } else { 
    ll <- cbind(as.matrix(x[, c("x_", "y_")]), 0)
  }
  
  ll <- ll * pi / 180
  xyz <- proj4::ptransform(ll, src.proj = p4, dst.proj = gproj)
  x$x_ <- xyz[,1]
  x$y_ <- xyz[,2]
  x$z_ <- xyz[,3]
  x
}



library(raster)

## Loading required package: sp

library(quadmesh)
library(geometry)

## Loading required package: magic

## Loading required package: abind

## 
## Attaching package: 'magic'

## The following object is masked from 'package:raster':
## 
##     shift

library(geosphere)
library(proj4)
library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:raster':
## 
##     intersect, select, union

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

mesh <- as.data.frame(geosphere::randomCoordinates(1024)) %>% 
  transmute(x_ = lon, y_ = lat, z_ = 0) %>% 
  globe(p4 = "+proj=longlat +datum=WGS84") 
index <- delaunayn(as.matrix(mesh))

## 
##      PLEASE NOTE:  As of version 0.3-5, no degenerate (zero area) 
##      regions are returned with the "Qt" option since the R 
##      code removes them from the triangulation. 
##      See help("delaunayn").

library(rgl)
rgl.triangles(as.matrix(mesh[t(surf.tri(as.matrix(mesh), delaunayn(as.matrix(mesh)))), ]), 
              col = rep(c('black', "grey"), each = 3))
rglwidget()