r - 计算并绘制任意 rasterLayer 的矢量场
问题描述
问题陈述:
我们可以绘制矢量场,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 单位)来计算和绘制矢量场。源代码在这里。
library(raster)
library(rasterVis)
r <- getData('alt', country='FRA', mask=TRUE)
r <- aggregate(r, 20)
vectorplot(r, par.settings=RdBuTheme())
结论:
回顾一下,我想取任意rasterLayer
高程,将其转换为 a 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))
解决方案
实际上,您要问的是将 2D scalar field转换为vector field。有几种不同的方法可以做到这一点。
raster 包包含函数terrain
,它创建新的栅格图层,该图层将为您提供所需矢量在每个点的角度(即aspect)及其大小(slope)。我们可以使用一点三角函数将它们转换为所使用的南北和东西基向量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://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
xend
yend
stat_quiver
xend
yend
geom_segment
u
v
推荐阅读
- r - R shiny-app DT:使用大于应用过滤器
- php - 蛋糕php中未定义的索引X-Api-Key
- osgi - 在新版本的 JBossFuse 中捆绑中突然出现未解决的约束
- api - 我如何自己构建 tensorflow-serving-api?
- assembly - 汇编 att x86-64 使用 C fopen 打开文件
- javascript - 无法对这个特定服务“未定义”目标命名空间进行 SOAP 调用
- java - Android 不知道如何修复这个索引越界异常 FirebaseRecyclerview
- c - 指针算术以 Code139 终止程序(分段错误)
- php - Twilio 异常在 codeigniter 中给出错误
- ios - 自动布局和底线文本字段