Tamponnage euclidien et géodésique en R

Pour comprendre la mise en mémoire tampon géodésique , l'équipe de développement du géotraitement d'Esri fait la distinction entre la mise en mémoire tampon euclidienne et géodésique. Ils concluent avec "la mise en mémoire tampon euclidienne effectuée sur les classes d'entités projetées peut produire des tampons trompeurs et techniquement incorrects. Cependant, la mise en mémoire tampon géodésique produira toujours des résultats qui sont géographiquement précis car les tampons géodésiques ne sont pas affectés par les distorsions introduites par les systèmes de coordonnées projetées".

Je dois travailler avec un jeu de données global ponctuel et les coordonnées ne sont pas projetées ( +proj=longlat +ellps=WGS84 +datum=WGS84). Existe-t-il une fonction pour créer un tampon géodésique dans R lorsque la largeur est donnée en unités métriques? Je suis au courant gBufferdu rgeospaquet. Cette fonction crée un tampon en unités de l'objet géographique utilisé ( exemple ), donc, je dois projeter les coordonnées pour pouvoir créer un tampon de X km souhaité. Projeter puis appliquer un gBuffermoyen fait en fait un tampon euclidien par opposition à un tampon géodésique dont j'ai besoin. Voici un code pour illustrer mes préoccupations:

require(rgeos)
require(sp)
require(plotKML)

# Generate a random grid-points for a (almost) global bounding box
b.box <- as(raster::extent(120, -120, -60, 60), "SpatialPolygons")
proj4string(b.box) <- "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"
set.seed(2017)
pts <- sp::spsample(b.box, n=100, type="regular")
plot(pts@coords)

# Project to Mollweide to be able to apply buffer with `gBuffer` 
# (one could use other projection)
pts.moll <- sp::spTransform(pts, CRSobj = "+proj=moll")
# create 1000 km buffers around the points
buf1000km.moll <- rgeos::gBuffer(spgeom = pts.moll, byid = TRUE, width = 10^6)
plot(buf1000km.moll)
# convert back to WGS84 unprojected
buf1000km.WGS84 <- sp::spTransform(buf1000km.moll, CRSobj = proj4string(pts))
plot(buf1000km.WGS84) # distorsions are present
# save as KML to better visualize distorted Euclidian buffers on Google Earth
plotKML::kml(buf1000km.WGS84, file.name = "buf1000km.WGS84.kml")

L'image ci-dessous représente les tampons euclidiens déformés (rayon de 1000 km) produits avec le code ci-dessus.

Robert J. Hijmans dans Introduction au paquet «géosphère» , la section 4 Point at distance and bearingdonne un exemple de la façon de faire des «polygones circulaires avec un rayon fixe, mais en coordonnées de longitude / latitude», qui je pense peut être appelé un «tampon géodésique». J'ai suivi cette idée et j'ai écrit du code qui, espérons-le, fait la bonne chose, mais je me demande s'il y a déjà une fonction R de tampon géodésique dans un paquet qui permet un rayon métrique en entrée:

require(geosphere)

make_GeodesicBuffer <- function(pts, width) {
    ### A) Construct buffers as points at given distance and bearing
    # a vector of bearings (fallows a circle)
    dg <- seq(from = 0, to = 360, by = 5)

    # Construct equidistant points defining circle shapes (the "buffer points")
    buff.XY <- geosphere::destPoint(p = pts, 
                                    b = rep(dg, each = length(pts)), 
                                    d = width)

    ### B) Make SpatialPolygons
    # group (split) "buffer points" by id
    buff.XY <- as.data.frame(buff.XY)
    id  <- rep(1:length(pts), times = length(dg))
    lst <- split(buff.XY, id)

    # Make SpatialPolygons out of the list of coordinates
    poly   <- lapply(lst, sp::Polygon, hole = FALSE)
    polys  <- lapply(list(poly), sp::Polygons, ID = NA)
    spolys <- sp::SpatialPolygons(Srl = polys, 
                                  proj4string = CRS(as.character("+proj=longlat +ellps=WGS84 +datum=WGS84")))
    # Disaggregate (split in unique polygons)
    spolys <- sp::disaggregate(spolys)
    return(spolys)
}

buf1000km.geodesic <- make_GeodesicBuffer(pts, width=10^6)
# save as KML to visualize geodesic buffers on Google Earth
plotKML::kml(buf1000km.geodesic, file.name = "buf1000km.geodesic.kml")

L'image ci-dessous représente les tampons géodésiques (1000 km de rayon).

Edit 2019-02-12 : Pour plus de commodité, j'ai encapsulé une version de la fonction dans le package geobuffer . N'hésitez pas à contribuer avec des demandes de tirage.

Dans la plupart des cas, il sera suffisamment précis pour utiliser un modèle sphérique de la terre - et le codage sera plus facile et les calculs beaucoup plus rapides.

Suivant les suggestions de M. Kennedy dans un commentaire, cette solution tamponne le pôle Nord (ce qui est facile: la limite du tampon se trouve à une latitude fixe), puis fait pivoter ce tampon dans n'importe quel emplacement souhaité.

La rotation est effectuée en convertissant le tampon d'origine en coordonnées géocentriques cartésiennes (XYZ), en faisant tourner celles avec une multiplication matricielle (rapide) le long du méridien principal vers la latitude cible, en reconvertissant ses coordonnées en géographiques (lat-lon), puis en les faisant tourner le tampon autour de l'axe de la Terre simplement en ajoutant la longitude cible à chaque seconde coordonnée.

Pourquoi le faire en deux étapes alors que (normalement) une seule multiplication matricielle fonctionnerait? Parce qu'il n'y a pas besoin de code spécial pour identifier les ruptures au méridien de +/- 180 degrés. Au lieu de cela, cette approche peut générer des longitudes au-delà de la plage d'origine (que ce soit -180 à 180 degrés ou 0 à 360 ou autre), mais ce faisant, les procédures standard de dessin de polygones (et d'autres procédures analytiques) fonctionneront correctement sans modification. Si vous devez avoir des longitudes dans une plage donnée, réduisez-les simplement à 360 degrés modulo à la toute fin: c'est rapide et facile.

Lors de la mise en mémoire tampon des points, tous les tampons ont généralement le même rayon. Cette solution modulaire permet une accélération dans ce cas: vous pouvez tamponner le pôle Nord puis le convertir en coordonnées XYZ une fois pour toutes. La mise en mémoire tampon de chaque point nécessite ainsi une multiplication matricielle (très rapide), une reconversion en coordonnées lat-lon et un décalage des longitudes (également très rapide). Attendez-vous à générer environ 10 000 tampons haute résolution (360 sommets) par seconde.

Ce Rcode montre les détails. Puisque son but est l'illustration, il n'utilise aucun module complémentaire; rien n'est caché ou enterré. Il comprend un test dans lequel un ensemble de points aléatoires est généré, mis en mémoire tampon et affiché en utilisant ses coordonnées lat-lon (géographiques) brutes. Voici la sortie:

degrees.to.radians <- function(phi) phi * (pi / 180)
radians.to.degrees <- function(phi) phi * (180 / pi)
#
# Create a 3X3 matrix to rotate the North Pole to latitude `phi`, longitude 0.
# Solution: A rotation is a linear map, and therefore is determined by its
#           effect on a basis.  This rotation does the following:
#           (0,0,1) -> (cos(phi), 0, sin(phi))  {North Pole (Z-axis)}
#           (0,1,0) -> (0,1,0)                  {Y-axis is fixed}
#           (1,0,0) -> (sin(phi), 0, -cos(phi)) {Destination of X-axis}
#
rotation.create <- function(phi, is.radians=FALSE) {
  if (!is.radians) phi <- degrees.to.radians(phi)
  cos.phi <- cos(phi)
  sin.phi <- sin(phi)
  matrix(c(sin.phi, 0, -cos.phi, 0, 1, 0, cos.phi, 0, sin.phi), 3)
}
#
# Convert between geocentric and spherical coordinates for a unit sphere.
# Assumes `latlon` in degrees.  It may be a 2-vector or a 2-row matrix.
# Returns an array with three rows for x,y,z.
#
latlon.to.xyz <- function(latlon) {
  latlon <- degrees.to.radians(latlon)
  latlon <- matrix(latlon, nrow=2)
  cos.phi <- cos(latlon[1,])
  sin.phi <- sin(latlon[1,])
  cos.lambda <- cos(latlon[2,])
  sin.lambda <- sin(latlon[2,])
  rbind(x = cos.phi * cos.lambda,
        y = cos.phi * sin.lambda,
        z = sin.phi)
}
xyz.to.latlon <- function(xyz) {
  xyz <- matrix(xyz, nrow=3) 
  radians.to.degrees(rbind(phi=atan2(xyz[3,], sqrt(xyz[1,]^2 + xyz[2,]^2)),
                           lambda=atan2(xyz[2,], xyz[1,])))
}
#
# Create a circle of radius `r` centered at the North Pole, oriented positively.
# `r` is measured relative to the sphere's radius `R`.  For the authalic Earth,
# r==1 corresponds to 6,371,007.2 meters.
#
# `resolution` is the number of vertices to use in a polygonal approximation.
# The first and last vertex will coincide.
#
circle.create <- function(r, resolution=360, R=6371007.2) {
  phi <- pi/2 - r / R                       # Constant latitude of the circle
  resolution <- max(1, ceiling(resolution)) # Assures a positive integer
  radians.to.degrees(rbind(phi=rep(phi, resolution+1),
                           lambda=seq(0, 2*pi, length.out = resolution+1)))
}
#
# Rotate around the y-axis, sending the North Pole to `phi`; then
# rotate around the new North Pole by `lambda`.
# Output is in geographic (spherical) coordinates, but input points may be
# in Earth-centered Cartesian or geographic.
# No effort is made to clamp longitudes to a 360 degree range.  This can 
# facilitate later computations.  Clamping is easily done afterwards if needed:
# reduce the longitude modulo 360 degrees.
#
rotate <- function(p, phi, lambda, is.geographic=FALSE) {
  if (is.geographic) p <- latlon.to.xyz(p)
  a <- rotation.create(phi)   # First rotation matrix
  q <- xyz.to.latlon(a %*% p) # Rotate the XYZ coordinates
  q + c(0, lambda)            # Second rotation
}
#------------------------------------------------------------------------------#
#
# Test.
#
n <- 50                  # Number of circles
radius <- 1.5e6          # Radii, in meters
resolution <- 360
set.seed(17)             # Makes this code reproducible

#-- Generate random points for testing.
centers <- rbind(phi=(rbeta(n, 2, 2) - 1/2) * 180,
                 lambda=runif(n, 0, 360))

system.time({
  #-- Buffer the North Pole and convert to XYZ once and for all.
  p.0 <- circle.create(radius, resolution=resolution) 
  p <- latlon.to.xyz(p.0)

  #-- Buffer the set of points (`centers`).
  circles <- apply(centers, 2, function(center) 
    rotate(p, center[1], center[2]))

  #-- Convert into an array indexed by (latlon, vertex, point id).
  circles <- array(circles, c(2, resolution+1, n))
})
#
# Display the buffers (if there are not too many).
#
if (n <= 1000) {
  #-- Create a background map area and graticule.
  xlim <- range(circles[2,,]) # Extent of all longitudes in the buffers
  plot(xlim, c(-90, 90), type="n", xlim=xlim, ylim=c(-90,90), asp=1,
       xlab="Longitude", ylab="Latitude",
       main=paste(n, "Random Disks of Radius", signif(radius/1e3, 3), "Km"),
       sub="Centers shown with gray dots")
  abline(v=seq(-360, 720, by=45), lty=1, col="#d0d0d0")
  abline(h=seq(-90, 90, by=30), lty=1, col="#d0d0d0")

  #-- Display the buffers themselves.
  colors <- terrain.colors(n, alpha=1/3) # Vary their colors
  invisible(sapply(1:n, function(i) {
    polygon(circles[2,,i], circles[1,,i], col=colors[i])
  }))

  #-- Show the original points (and, optionally, labels).
  points(centers[2,], centers[1,], pch=21, bg="Gray", cex=min(1, sqrt(25/n)))
  # text(centers[2,], centers[1,], labels=1:n, cex=min(1, sqrt(100/n)))
}