Al comprender el almacenamiento en búfer geodésico , el equipo de desarrollo de geoprocesamiento de Esri distingue entre el almacenamiento en búfer geodésico y euclidiano. Concluyen con "El almacenamiento en búfer euclidiano realizado en clases de entidad proyectadas puede producir amortiguadores engañosos y técnicamente incorrectos. Sin embargo, el almacenamiento en búfer geodésico siempre producirá resultados geográficamente precisos porque los amortiguadores geodésicos no se ven afectados por las distorsiones introducidas por los sistemas de coordenadas proyectados".
Tengo que trabajar con un conjunto de datos global de puntos y las coordenadas no están proyectadas ( +proj=longlat +ellps=WGS84 +datum=WGS84
). ¿Existe una función para crear un búfer geodésico en R cuando el ancho se da en unidades métricas? Soy consciente de gBuffer
de rgeos
paquete. Esta función crea un búfer en unidades del objeto espacial que se usa ( ejemplo ), por lo tanto, tengo que proyectar las coordenadas para poder crear un búfer de X km deseado. Proyectando y luego aplicando un gBuffer
medio en realidad haciendo un búfer euclidiano en lugar de uno geodésico que necesito. A continuación hay un código para ilustrar mis preocupaciones:
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")
La imagen a continuación muestra los buffers euclidianos distorsionados (radio de 1000 km) producidos con el código desde arriba.
Robert J. Hijmans en Introducción al paquete de "geosfera" , sección 4 Point at distance and bearing
da un ejemplo de cómo hacer "polígonos circulares con un radio fijo, pero en coordenadas de longitud / latitud", que creo que se puede llamar un "amortiguador geodésico". Seguí esta idea y escribí un código que con suerte hace lo correcto, pero me pregunto si ya hay alguna función R de búfer geodésico en algún paquete que permita el radio métrico como entrada:
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")
La siguiente imagen muestra los amortiguadores geodésicos (radio de 1000 km).
Editar 2019-02-12 : por conveniencia, incluí una versión de la función en el paquete de geobuffer . Siéntase libre de contribuir con solicitudes de extracción.
R
puedo hacer eso, es una gran sugerencia. Pero dado que para un modelo esférico de la Tierra esta es una proyección tan simple, es lo suficientemente simple como para escribir el código directamente.Respuestas:
Para la mayoría de los propósitos, será lo suficientemente preciso como para usar un modelo esférico de la Tierra, y la codificación será más fácil y los cálculos mucho más rápidos.
Siguiendo las sugerencias de M. Kennedy en un comentario, esta solución amortigua el Polo Norte (lo cual es fácil: el límite del buffer se encuentra en una latitud fija) y luego gira este buffer en cualquier ubicación deseada.
La rotación se efectúa convirtiendo el búfer original en coordenadas cartesianas geocéntricas (XYZ), rotando aquellas con una multiplicación de matriz (rápida) a lo largo del meridiano principal a la latitud objetivo, convirtiendo sus coordenadas nuevamente en Geográfico (lat-lon) y luego girando el búfer alrededor del eje de la Tierra simplemente agregando la longitud objetivo a cada segunda coordenada.
¿Por qué hacerlo en dos pasos cuando (normalmente) una sola multiplicación de matriz funcionaría? Porque no hay necesidad de un código especial para identificar los descansos en el meridiano de +/- 180 grados. En cambio, este enfoque puede generar longitudes más allá del rango original (ya sea -180 a 180 grados o 0 a 360 o lo que sea), pero al hacerlo, los procedimientos estándar de dibujo de polígonos (y otros procedimientos analíticos) funcionarán bien sin modificación. Si debe tener longitudes en un rango determinado, simplemente reduzca el módulo 360 grados al final: eso es rápido y fácil.
Cuando se almacenan puntos en el búfer, generalmente todos los búferes tienen el mismo radio. Esta solución modular permite una aceleración en este caso: puede amortiguar el Polo Norte y luego convertirlo en coordenadas XYZ de una vez por todas. Para amortiguar cada punto, se requiere una multiplicación de matriz (muy rápida), conversión de nuevo a coordenadas lat-lon y desplazamiento de las longitudes (también muy rápido). Espere generar alrededor de 10,000 buffers de alta resolución (360 vértices) por segundo.
Este
R
código muestra los detalles. Como su propósito es ilustrativo, no utiliza paquetes adicionales; Nada está escondido o enterrado. Incluye una prueba en la que se genera, almacena y almacena un conjunto de puntos aleatorios utilizando sus coordenadas lat-lon (geográficas) sin procesar. Aquí está la salida:fuente