计算和绘制任意栅格层的矢量场

Question

问题陈述：

使用ggquiver::geom_quiver（），我们可以绘制向量场，前提是我们知道x，y，xend和yend的值。

1. 我如何为任意高程的RasterLayer计算这些参数？
2. 我如何确保这些箭头的大小显示该向量的坡度，以便箭头在该位置的梯度成比例地不同长度显示（例如下面的第一张图）？

---

背景：

# 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单位）来计算和绘制矢量场。 源代码在这里。

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))

- Rich Pauloo

---

为什么是任意的，难道大小斜率不是按照栅格单元大小定向在方位方向上进行归一化的吗？然后从每个单元格的中心计算距离和方位角。请参见geosphere pkg。 - wildintellect

我认为这个答案会对你有所帮助：https://dev59.com/EnHYa4cB1Zd3GeqPOral#16371769 - Oscar Perpiñán

Answer 1

您实际上是要将一个二维标量场转换为向量场。有几种不同的方法可以做到这一点。

raster包中包含函数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)

然而，在大多数情况下，这不会产生一个好的图形。如果您不先聚合栅格，您将得到图像上每个像素的一个梯度，这样绘制效果不佳。另一方面，如果您聚合了，您将得到一个漂亮的矢量场，但您的栅格将看起来“块状”。因此，拥有单个数据框可能不是最好的选择。

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

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://istack.dev59.com/tXUXO.webp")

# 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值的大小。

- Allan Cameron