我想使用 R 中的 leaflet 包绘制一些空间数据,但是生成的光栅图像与参考网格相比似乎发生了偏移。我怀疑地图投影问题,但我不是该主题的专家,因此任何帮助将不胜感激。
这是绘制地图的最小代码:
library(leaflet)
library(sp)
library(raster)
set.seed(111)
# create dummy data -rectangular grid with random values
m = 10
n = 10
x = seq(45,48,length.out = m)
y = seq(15,18,length.out = n)
X = matrix(rep(x, each = n), nrow = n)
Y = matrix(rep(y, m), nrow = n)
# collector dataframe
points = data.frame(value = rnorm(n*m), lng = c(Y), lat = c(X))
## create raster grid
s = SpatialPixelsDataFrame(points[,c('lng', 'lat')], data = points)
# set WGS84 projection
crs(s) = sp::CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
r = raster(s)
# add coloring
pal = colorNumeric(c("#0C2C84", "#f7f7f7", "#F98009"), points$value,
na.color = "transparent")
## plot map
leaflet() %>% addProviderTiles("CartoDB.Positron") %>%
addRasterImage(r, colors = pal, opacity = 0.6)
This produces this map, which is ok at first sight:
如果在此地图上添加网格:
## grid
dx = diff(x)[1]/2
dy = diff(y)[1]/2
rect_lng = sapply(points$lng, function(t) c(t-dx, t+dx))
rect_lat = sapply(points$lat, function(t) c(t-dy, t+dy))
leaflet() %>% addProviderTiles("CartoDB.Positron") %>%
addRectangles(
lng1=rect_lng[1,], lat1=rect_lat[1,],
lng2=rect_lng[2,], lat2=rect_lat[2,],
fillColor = "transparent",
weight = 1
) %>%
addRasterImage(r, colors = pal, opacity = 0.6)
地图看起来像这样:
Here we can see that the grids are not matching.
这种不匹配的原因是什么?怎么可能消除呢?我尝试了各种预测,但没有成功。唯一有效的就是使用addRectangle
代替addRasterImage
,但这需要更多的计算并减慢进程,因此我想避免。请注意,在上面的示例中addRectangle
仅用于参考,在最终代码中我不想使用它。
对于具有更多单元格(网格)的地图,不匹配非常大,可能大于单个单元格的大小。
预先感谢您的任何帮助。
EDIT
该问题可能与椭球体和球体投影之间的投影问题有关,请参阅最后一个问题here http://proj.maptools.org/faq.html:
要在球体上的 WGS84 和墨卡托之间进行转换,将会有
Y 墨卡托坐标发生重大变化。这是因为
在内部,cs2cs 必须调整纬度/经度坐标
位于球体上到位于 WGS84 基准上,该基准具有相当大的
不同形状的椭球体。
但是,我无法使用推荐的“技巧”解决问题:+nadgrids=@null
.