计算和绘制任意栅格图层的向量场。

Posted

tags:

篇首语:本文由小常识网(cha138.com)小编为大家整理,主要介绍了计算和绘制任意栅格图层的向量场。相关的知识,希望对你有一定的参考价值。

问题说明:

随着 ggquiver::geom_quiver() 我们可以绘制向量场,只要我们知道 x, y, xendyend.

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

enter image description here

一个相关的应用程序使用 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())

enter image description here


结论。

回顾一下,我想以一种任意的方式 rasterLayer 的海拔高度,将其转换为 data.frame,计算出 x, y, xmaxymax 高程矢量场的分量,这些分量的箭头大小可以显示出该点的相对坡度(如上面的图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))
答案

实际上,你所要求的是将一个二维的 标量场 变成 向量场. 有几种不同的方法可以做到这一点。

栅格包中包含的功能是 terrain,创建新的光栅层,使您在每个点上都能得到所需矢量的角度(即 方面),以及其幅度(the 坡度). 我们可以用一点三角法将其转换为南北和东西基向量,用于 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)

enter image description here

现在要显示它在一个 任意 栅格(前提是它已经分配了一个投影),让我们在这幅转换为栅格的图像上测试一下。这一次,我们有一个较低的分辨率,所以我们选择一个较小的集合值。我们还将为最低值选择一个透明的颜色,以提供一个更好的图形。

enter image description here

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

enter image description here

* 我应该指出的是 geom_quiver的映射美学,将参数称为 uv,表示指向北和东的基向量。的 ggquiver 包将它们转换为 xendyend 值,使用 stat_quiver. 如果您喜欢使用 xendyend 值,你可以直接使用 geom_segment 来绘制你的向量场,但这使得控制箭头的外观更加复杂。因此,这个解决方案将找到箭头的大小。uv 值代替。

以上是关于计算和绘制任意栅格图层的向量场。的主要内容,如果未能解决你的问题,请参考以下文章

ArcGIS server10.1如何动态添加栅格图层

如何使用arcgis将栅格图转为矢量图,

混合和过度绘制(图层性能 15.3)

在Arcmap中如何创建地图图层

arcgis中如何计算dem图层单元栅格的曲面面积、和投影到平面的面积、就是要计算地表粗糙度

给定一个任意单位向量,计算任意正交单位向量的最佳方法是啥?