操作地图切片数据以与 R 中的栅格数据对齐

Question

我正在尝试使用 rgl 在 3D 表面上绘制地图图块，但无法弄清楚如何正确对齐数据。这可能与 R 在矩阵和栅格之间转换时添加 90 度旋转的行为有关，但我也发现需要在代码中添加翻转以获得正确的结果。工作流程有点棘手，所以我把它做成了一个函数来显示结果与输入变量的关系。演示：

require(raster)
require(akima)
require(OpenStreetMap)
require(rgl)

wgs84 = '+proj=longlat +datum=WGS84'

plot_3d_tile = function(z, xlims, ylims, zscale=1, zoom, crs, plot_rasters=F, ...){
  # specify raster's spatial info
  extent(z) = c(xlims, ylims); crs(z) = crs

  # extend range slightly to crop back to rect after reproj
  osm_x = extendrange(r=xlims); osm_y = extendrange(r=ylims)

  # get OSM map tile & reproject to wgs84
  m = raster(openproj(openmap(c(osm_y[2],osm_x[1]), c(osm_y[1],osm_x[2]), zoom=zoom)))
  m = crop(flip(m,'y'), extent(z))             # FLIPPED
  if(plot_rasters) plotRGB(m)

  # coerce to lists of points for akima::interp
  pts_m = rasterToPoints(m); pts_z = rasterToPoints(z)

  # resizes z to match tile
  intp = interp(x=pts_z[,1], y=pts_z[,2], z=pts_z[,3], 
                xo=unique(pts_m[,1]), yo=unique(pts_m[,2]))

  # get matrix of interpolated z values and convert back to spatial
  z2 = flip(raster(apply(intp$z, 1, rev)),'y') # FLIPPED AND ROTATED
  cat("dimensions match? ", dim(z2) == dim(raster(m)))    # check dimensions match up
  extent(z2) = extent(xlims, ylims); crs(z2) = crs  # spatialise
  if(plot_rasters) plot(z2, asp=T)
  pts_z2 = rasterToPoints(z2)

  # create hex colour vector from tile values
  col_data = getValues(m)
  cols = rgb(col_data[,1], col_data[,2], col_data[,3], maxColorValue = 255)

  # plot 3d extruded map tile
  rgl.open(); bg3d("white")
  rgl.surface(unique(pts_z2[,1]), unique(pts_z2[,2]), pts_z2[,3]*zscale, 
              color=cols, specular="black",  back="lines", asp=T, ...)
  results <<- list(z=z, z2=z2, m=m, intp=intp, pts_m=pts_m, pts_z=pts_z, pts_z2=pts_z2)
}

z1 = raster(volcano)
xlims = c(-0.24, -0.1)
ylims = c(51.4, 51.58)

现在来测试一下：

plot_3d_tile(z1, xlims, ylims, zscale=1/3000, zoom=9, crs=wgs84)

plot_3d_tile(z1, xlims, ylims, zscale=1/3000, zoom=10, crs=wgs84)

plot_3d_tile(z1, xlims, ylims, zscale=1/3000, zoom=11, crs=wgs84)

如您所见，随着 OSM 缩放级别的增加，开始看起来不错的东西会逐渐变形。我怀疑就翻转和旋转而言我有问题，但这是迄今为止我实现的最接近的组合。我知道一个特定的问题，但代码可能对其他人有用，所以我在这里发布。提前致谢。

Answer 1

我想出了一个解决方法。除了将颜色信息输入rgl.surface 之外，还可以为其texture 参数提供一个png。这可能意味着我使表面栅格与图块尺寸匹配的代码有点多余，尽管如果表面分辨率低得多，它仍然可以实现更平滑的插值。工作功能：

plot_3d_tile = function(z, xlims, ylims, zscale=1, zoom, crs, plot_rasters=F, ...){
  # specify raster's spatial info
  extent(z) = c(xlims, ylims); crs(z) = crs
  if(plot_rasters) plot(z, asp=T, main='z')

  # extend range slightly to crop back to rect after reproj
  osm_x = extendrange(r=xlims); osm_y = extendrange(r=ylims)

  # get OSM map tile & reproject to wgs84
  m = raster(openproj(openmap(c(osm_y[2],osm_x[1]), c(osm_y[1],osm_x[2]), zoom=zoom)))
  m = crop(flip(m,'y'), extent(z))             # FLIPPED
  if(plot_rasters) plotRGB(m, asp=T, main='m')
  png('plot.png', width=ncol(m), height=nrow(m))
  plotRGB(m)
  dev.off()

  # coerce to lists of points for akima::interp
  pts_m = rasterToPoints(m); pts_z = rasterToPoints(z)

  # resizes z to match tile
  intp = interp(x=pts_z[,1], y=pts_z[,2], z=pts_z[,3], 
                xo=unique(pts_m[,1]), yo=unique(pts_m[,2]))

  # get matrix of interpolated z values and convert back to spatial
  z2 = raster(apply(intp$z, 1, rev)) # ROTATED
  cat("dimensions match? ", dim(z2) == dim(raster(m)))    # check dimensions match up
  extent(z2) = extent(xlims, ylims); crs(z2) = crs  # spatialise
  if(plot_rasters) plot(z2, asp=T, main='z2')
  pts_z2 = rasterToPoints(z2)

  # plot 3d extruded map tile
  bg3d("white")
  rgl.surface(unique(pts_z2[,1]), unique(pts_z2[,2]), pts_z2[,3]*zscale, 
              texture='plot.png', specular="black",  back="lines", asp=T, ...)
  results <<- list(z=z, z2=z2, m=m, intp=intp, pts_m=pts_m, pts_z=pts_z, pts_z2=pts_z2)
}

rgl.open()

plot_3d_tile(z1, xlims, ylims, zscale=1/3000, zoom=12, crs=wgs84, plot_rasters=T)