计算和绘制任意 rasterLayer 的矢量场

Question

问题陈述：

使用 ggquiver::geom_quiver() 我们可以绘制向量场，前提是我们知道 x、y、xend 和 yend。

如何计算任意 RasterLayer 海拔的这些参数？
我如何确保这些箭头的大小指示特定向量的斜率，以便箭头显示不同的长度与该位置的梯度成比例（例如，下面的第一个图）？

背景：

# ggquiver example

library(tidyverse)
library(ggquiver)
expand.grid(x=seq(0,pi,pi/12), y=seq(0,pi,pi/12)) %>%
  ggplot(aes(x=x,y=y,u=cos(x),v=sin(y))) +
  geom_quiver()

相关方法使用 rasterVis::vectorplot，它依赖于 raster::terrain（假设场单位 == CRS 单位）来计算和绘制矢量场。 Source code is here.

library(raster)
library(rasterVis)
r <- getData('alt', country='FRA', mask=TRUE)
r <- aggregate(r, 20)
vectorplot(r, par.settings=RdBuTheme())

结论：

为了复习，我想获取任意 rasterLayer 高程，将其转换为 data.frame，计算 x、y、xmax 和 ymax 高程矢量场的分量，这些分量调整箭头的大小以显示该点的相对坡度（如上图 1 和 2 所示），并使用 ggquiver 绘制。类似于：

names(r) <- "z"
rd <- as.data.frame(r, xy=TRUE)

# calculate x, y, xend, yend for gradient vectors, add to rd, then plot

ggplot(rd) + 
  geom_raster(aes(x, y, fill = z)) + 
  geom_quiver(aes(x, y, xend, yend))

Answer 1

实际上您要求的是转换 2D scalar field into a vector field。有几种不同的方法可以做到这一点。

光栅包包含函数 terrain，它会创建新的光栅图层，为您提供每个点所需矢量的角度（即 纵横比 )，及其大小（斜率）。我们可以使用一点三角函数将它们转换为 ggquiver 使用的南北向和东西向基向量，并将它们添加到我们的原始栅格中，然后再将整个东西变成数据框。*

terrain_raster <- terrain(r, opt = c('slope', 'aspect'))
r$u <- terrain_raster$slope[] * sin(terrain_raster$aspect[])
r$v <- terr$slope[] * cos(terr$aspect[])
rd <- as.data.frame(r, xy = TRUE)

然而，在大多数情况下，这不会成为一个好的情节。如果您不首先聚合栅格，图像上的每个像素都会有一个梯度，这将无法很好地绘制。另一方面，如果您做聚合，您将有一个很好的矢量场，但您的栅格看起来 "blocky"。因此，为您的绘图使用单个数据框可能不是最好的方法。

以下函数将获取栅格并使用叠加矢量场绘制它。您可以在不影响栅格的情况下调整矢量场的聚合程度，并且可以为栅格指定任意颜色矢量。

raster2quiver <- function(rast, aggregate = 50, colours = terrain.colors(6))
{
  names(rast) <- "z"
  quiv <- aggregate(rast, aggregate)
  terr <- terrain(quiv, opt = c('slope', 'aspect'))
  quiv$u <- terr$slope[] * sin(terr$aspect[])
  quiv$v <- terr$slope[] * cos(terr$aspect[])
  quiv_df <- as.data.frame(quiv, xy = TRUE)
  rast_df <- as.data.frame(rast, xy = TRUE)

  print(ggplot(mapping = aes(x = x, y = y, fill = z)) + 
          geom_raster(data = rast_df, na.rm = TRUE) + 
          geom_quiver(data = quiv_df, aes(u = u, v = v), vecsize = 1.5) +
          scale_fill_gradientn(colours = colours, na.value = "transparent") +
          theme_bw())

  return(quiv_df)
}

因此，在您的法国示例中尝试一下，在首先定义相似的调色板之后，我们得到

pal <- c("#B2182B", "#E68469", "#D9E9F1", "#ACD2E5", "#539DC8", "#3C8ABE", "#2E78B5")

raster2quiver(getData('alt', country = 'FRA', mask = TRUE), colours = pal)

现在为了证明它适用于任意光栅（前提是它已指定投影）让我们在转换为光栅的图像上对其进行测试。这一次，我们有一个较低的分辨率，所以我们选择一个较小的聚合值。我们还将为最低值选择透明颜色以提供更好的绘图：

rast <- raster::raster("https://i.stack.imgur.com/tXUXO.png")

# Add a fake arbitrary projection otherwise "terrain()" doesn't work:
projection(rast) <- "+proj=lcc +lat_1=48 +lat_2=33 +lon_0=-100 +ellps=WGS84"

raster2quiver(rast, aggregate = 20, colours = c("#FFFFFF00", "red"))

*_{我应该指出，geom_quiver 的映射美学采用名为 u 和 v 的参数，它们表示基向量指向北方和东方。 ggquiver 包使用 stat_quiver 将它们转换为 xend 和 yend 值。如果您更喜欢使用 xend 和 yend 值，您可以只使用 geom_segment 来绘制矢量场，但这会使控制箭头的外观变得更加复杂。因此，此解决方案将找到 u 和 v 值的大小。}

计算和绘制任意 rasterLayer 的矢量场

calculate and plot vector field of an arbitrary rasterLayer

r

geospatial

ggplot2

r-raster

问题陈述：

背景：

结论：