栅格和多边形未对齐，即使范围和 projargs 相同答案

【问题标题】：Raster and polygon not aligning, even when extent and projargs are same栅格和多边形未对齐，即使范围和 projargs 相同
【发布时间】：2019-02-22 12:12:19
【问题描述】：

我正在尝试将空间多边形（源自 icosa 生成的网格）覆盖在全局栅格地图上。我已确保对象的边界框/范围相等，并且它们都投影到相同的 CRS，但多边形与栅格的范围不正确对齐。有任何想法吗？（此处提供基本图像：https://imgur.com/a/c1oVm6y）

# required packages
library(icosa)
library(raster)
library(rgdal)
library(rgeos)

# define projection string
equi <- CRS("+proj=longlat +datum=WGS84")

# load data and set extent
lad <- raster("Data/Maps/Ladinian_grey.jpg", crs = equi)
extent(lad) = c(-180, 180, -90, 90)

# generate grid and convert to spatial polygon
large_grid <- hexagrid(c(2, 3))
large_grid <- newsp(large_grid)
map_grid <- SpPolygons(large_grid, res = 50)
map_grid <- spTransform(map_grid, equi)

# plot
plot(lad)
plot(map_grid, col = NA, add = TRUE)

【问题讨论】：

请提供一个可重现的例子。 equi 是什么？您在raster( , crs=equi) 中使用它可能是错误的。 newsp 是什么？你为什么这样做：map_grid@bbox[1,1] <- extent(lad)[1] 等？那是不对的;永远不要手动更改这些值。
我忘记在我的原始代码中定义“equi”（它被用作 equirectangular proj4 参数的存储）。我在示例中对此进行了更新。 'newsp' 是 icosa 包中的一个函数。为了清晰起见，新代码已更新，尽管它的功能与我的原始代码相同。至于 bbox，我假设空间多边形的边界框和栅格的范围是类似的，因此尝试用后者更新前者的边界，以便它们的尺寸匹配。这样做并没有解决我的问题，但我按照我的匹配对象尺寸的假设保留了它
您不应更改 bbox 值，因为它们描述了对象的范围，但无论如何都不要更改它。所以当你这样做时你会弄得一团糟。相反，您可以查看值是什么，例如：extent(map_grid)
另外，这里可能无关紧要，但"+proj=longlat +datum=WGS84" 不是“等直角”。那将是proj=equi。
我进一步编辑了我的代码以反映以下内容。我看了看，范围确实不同，但是没有用于定义空间多边形范围的正式插槽。有没有其他方法可以编辑其尺寸以使其与光栅的尺寸对齐。无论字符串是什么，它们都被投影到同一个 CRS。最终我将两者都转换为 robinson CRS，但对齐问题仍然存在，所以为了简单起见，我将 robinson 投影排除在外。

标签： r polygon

【解决方案1】：

很难看出没有数据会发生什么，以及你所做的奇怪事情会改变对象的范围。下面的示例表明它应该可以工作：

library(raster)
library(icosa)

r <- raster()
values(r) <- 1:ncell(r)
large_grid <- hexagrid(c(2, 3))
large_grid <- newsp(large_grid)
map_grid <- SpPolygons(large_grid, res = 50)
plot(r)
lines(map_grid)

您显示的错误对齐可能是由绘制后调整绘图大小引起的。

它也适用于您的数据和（简化的）代码：

library(icosa)
library(raster)

lad <- raster("f91xQIt.jpg")
extent(lad) = c(-180, 180, -90, 90)
crs(lad) <- "+proj=longlat +datum=WGS84"

large_grid <- hexagrid(c(2, 3))
large_grid <- newsp(large_grid)
map_grid <- SpPolygons(large_grid, res=50)

# I use image so that it is clearer where all the lines come from
image(lad)
plot(map_grid, col = NA, add = TRUE)

而且程度符合预期

map_grid
#class       : SpatialPolygons 
#features    : 362 
#extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#crs         : +proj=longlat +a=6371007 +b=6371007

【讨论】：