在 R 中对齐栅格和 shapefile 的范围：从具有相同 CRS 的 QGIS 导入

Question

我有一个小光栅。它是来自 Digital Globe 底图的单个图块（缩放 = 16）；它是标准的 256 x 256 像素。在 QGIS 中，我创建了一个多边形 shapefile 并添加了大约 50 个特征。栅格和 shapefile 在同一个 CRS 中。在 QGIS 中，它们可以很好地对齐，如下图所示。

但是，当我在 R 中同时打开光栅和 shapefile 时，两者不再对齐。同样，两者都在同一个 CRS 中。可重现代码如下，shapefile 和光栅可以从 GitHub 文件夹here 下载。我正在使用 raster 包中的 brick() 来保持 RGB 波段分开。

#load
library(raster)
img <- brick("26261.png")
proj4string(img) <- CRS("+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs")
trainData <- shapefile("26261_train.shp")

#shape plots fine alone
plot(trainData, axes = TRUE)
#raster plots fine alone
plot(img,1)
#two do not match
plot(trainData, add = TRUE)

#summary
img
trainData 

问题似乎是由于两个文件具有不同的范围，如摘要所示：

> img
class       : RasterBrick 
dimensions  : 256, 256, 65536, 4  (nrow, ncol, ncell, nlayers)
resolution  : 1, 1  (x, y)
extent      : 0, 256, 0, 256  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs 
data source : /.../26261.png 
names       : X26261.1, X26261.2, X26261.3, X26261.4 
min values  :        0,        0,        0,        0 
max values  :      255,      255,      255,      255 

> trainData
class       : SpatialPolygonsDataFrame 
features    : 49 
extent      : 4.38, 244.9, -232.72, -6.48  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs 
variables   : 3
names       : id, cat1, cat2 
min values  : NA,    1,    0 
max values  : NA,    7,    2 

但是，当我在 QGIS 中检查两个文件的范围时，它们重叠。下图是 QGIS 中的栅格。 shapefile 和栅格的 y 范围都在 [-256, 0] 范围内。

brick() 函数似乎将 x 和 y 范围都设置为正值，尽管我已经查看了包/函数文档并且不明白为什么会这样。如何让这两个文件在 R 中对齐？

我尝试将栅格导出为 GeoTIFF，但这也不起作用。我还尝试使用其他一些 R 包（如 rgdal）读取 shapefile。如果在 QGIS 中导出 shapefile，以栅格为中心并选择“地图视图范围”，我可以让它“工作”，但这不是最佳解决方案，并且（a）我需要使用地图图块并且不想手动缩放到每个图块，并且（b）这并不能解释我所犯的错误。我什至不完全确定我的问题是从 QGIS 导出还是导入到 R，尽管我很确定我的错误出在 brick() 中。

注意：有一个类似的问题here，但尽管有问题的标题，但错误在于坐标参考系。

Answer 1

感谢@LoBu 的评论，虽然没有正确地进行地理参考并不能解释为什么从 QGIS 到 R 的导出会通过 brick() 翻转栅格的 y 范围，但在 R 中正确地对地图图块进行地理参考确实可以解决更大的问题。

这些图块是使用 QGIS 的 Qtiles 插件从 Digital Globe 底图中获取的“滑动地图”图块。因此，它们的索引如here 所述，对于 TMS 磁贴，实现仅略有不同。 Slippy 地图瓦片使用 NW 角作为参考，而 TMS 瓦片使用 SW 角作为参考。我们可以使用目录结构和文件名来地理参考每个图块。

对于我的问题中给出的示例图块，以下是 R 中的一个实现，改编了上一个链接中的 Python 示例：

#funtion to take tile names and return extent
num2deg <- function(xtile,ytile,zoom){

  n = 2^zoom

  lon_deg_nw = xtile / n * 360.0 - 180.0
  lat_rad_nw = atan(sinh(pi * (1 - 2 * ytile / n)))
  lat_deg_nw = lat_rad_nw * 180 / pi

  lon_deg_se = (xtile + 1) / n * 360.0 - 180.0
  lat_rad_se = atan(sinh(pi * (1 - 2 * (ytile - 1) / n)))
  lat_deg_se = lat_rad_se * 180 / pi

  #returns extent: vector (length=4; order= xmin, xmax, ymin, ymax)
  return(c(lon_deg_nw, lon_deg_se, lat_deg_se, lat_deg_nw))
}

#example get tile 
library(raster)
tile <- brick(".../16/39281/26261.png")

#get extent for raster
my.extent <- num2deg(xtile = 39281, ytile = 26261, zoom = 16)

#georeference and plot  
extent(tile) <- my.extent
projection(tile) <- CRS("+proj=longlat +datum=WGS84")
plot(tile, 1) 

仅绘制栅格的第一个波段，并对其进行了正确的地理参考：

我们还可以索引 Qtiles 导出的目录和子目录，以便我们可以根据需要检索和地理参考一个或多个图块：

library(stringr)
library(reshape2)

# index directories of slippy map tiles
# tiles obtain with Qtiles plugin for zoom levels 10-16
x <- list.files(".../tiles/", recursive = TRUE) 
index <- strsplit(x,"/") %>%
  data.frame() %>% 
  t() %>% 
  data.frame() 
# get only zoom level 16
index <- subset(index, X1 == "16")
# x and y array of tile reference numbers
index <- dcast(index, X2 ~ X3, value.var = "X3")
x.index <- as.character(index[,1])
y.index <- as.character(index[1,2:(length(index[1,]))]) 
y.index <- gsub(".png","",y.index)
x.index <- as.numeric(x.index)
y.index <- as.numeric(y.index)

#example usage
my.extent <- num2deg(xtile = x.index[1], ytile = y.index[1], zoom = 16)