Em Entendendo o buffer geodésico , a equipe de desenvolvimento de geoprocessamento da Esri distingue entre buffer euclidiano e geodésico. Eles concluem com "O buffer euclidiano realizado em classes de recursos projetados pode produzir buffers enganosos e tecnicamente incorretos. No entanto, o buffer geodésico sempre produzirá resultados geograficamente precisos porque os buffers geodésicos não são afetados pelas distorções introduzidas pelos sistemas de coordenadas projetadas".
Eu tenho que trabalhar com um conjunto de dados global pontual e as coordenadas não são projetadas ( +proj=longlat +ellps=WGS84 +datum=WGS84
). Existe uma função para criar um buffer geodésico em R quando a largura é fornecida em unidades métricas? Estou ciente gBuffer
do rgeos
pacote. Esta função cria um buffer em unidades do objeto espacial que é usado ( exemplo ), então, eu tenho que projetar as coordenadas para poder criar um buffer de X km desejado. Projetar e, em seguida, aplicar um gBuffer
meio, na verdade, fazer um buffer euclidiano em vez de um geodésico que eu preciso. Abaixo está um código para ilustrar minhas preocupações:
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")
A imagem abaixo mostra os buffers euclidianos distorcidos (raio de 1000 km) produzidos com o código acima.
Robert J. Hijmans, na introdução ao pacote "geosfera" , seção 4 Point at distance and bearing
fornece um exemplo de como criar "polígonos circulares com um raio fixo, mas em coordenadas de longitude / latitude", que eu acho que podem ser chamadas de "buffer geodésico". Adotei essa ideia e escrevi um código que, com sorte, faz a coisa certa, mas me pergunto se já existe alguma função R do buffer geodésico em algum pacote que permita o raio 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")
A imagem abaixo mostra os buffers geodésicos (raio de 1000 km).
Edit 2019-02-12 : Por conveniência, envolvi uma versão da função no pacote geobuffer . Sinta-se livre para contribuir com solicitações de recebimento.
R
pode fazer isso - é uma ótima sugestão. Mas como para um modelo esférico da Terra essa é uma projeção tão simples, é simples o suficiente escrever o código diretamente.Respostas:
Para a maioria dos propósitos, será preciso o suficiente para usar um modelo esférico da Terra - e a codificação será mais fácil e os cálculos muito mais rápidos.
Seguindo as sugestões de M. Kennedy em um comentário, esta solução armazena em buffer o Pólo Norte (o que é fácil: o limite do buffer fica em uma latitude fixa) e depois gira esse buffer para qualquer local desejado.
A rotação é efetuada convertendo o buffer original em coordenadas cartesianas geocêntricas (XYZ), girando aquelas com uma multiplicação de matriz (rápida) ao longo do Meridiano Principal para a latitude alvo, convertendo suas coordenadas de volta para Geográfica (lat-lon) e girando o buffer ao redor do eixo da Terra simplesmente adicionando a longitude alvo a cada segunda coordenada.
Por que fazê-lo em duas etapas quando (normalmente) uma multiplicação de matriz única funcionaria? Porque não há necessidade de código especial para identificar as quebras no meridiano de +/- 180 graus. Em vez disso, essa abordagem pode gerar longitudes além do intervalo original (-180 a 180 graus ou 0 a 360 ou o que for), mas ao fazer isso, os procedimentos padrão de desenho de polígonos (e outros procedimentos analíticos) funcionarão bem sem modificação. Se você precisa ter longitudes em um determinado intervalo, basta reduzi-las em 360 graus no final: isso é rápido e fácil.
Ao armazenar pontos em buffer, normalmente todos os buffers têm o mesmo raio. Esta solução modular permite uma aceleração neste caso: você pode armazenar em buffer o Polo Norte e convertê-lo em coordenadas XYZ de uma vez por todas. O buffer de cada ponto exige, assim, a multiplicação (muito rápida) da matriz, a conversão de volta para as coordenadas lat-lon e a mudança das longitudes (também muito rápidas). Espere gerar cerca de 10.000 buffers de alta resolução (360 vértices) por segundo.
Este
R
código mostra os detalhes. Como seu objetivo é ilustrativo, ele não usa pacotes complementares; nada está oculto ou enterrado. Ele inclui um teste no qual um conjunto de pontos aleatórios é gerado, armazenado em buffer e exibido usando suas coordenadas lat-lon (Geographic) brutas. Aqui está a saída:fonte