计算和绘制任意栅格图层的矢量场 [英] calculate and plot vector field of an arbitrary rasterLayer

问题描述

问题陈述:

使用 ggquiver :: geom_quiver()，我们可以绘制矢量场，只要我们知道 x ， y ， xend 和 yend .

1. 如何为任意 RasterLayer 高程计算这些参数?
2. 如何确保这些箭头的大小指示特定矢量的斜率，以使箭头在该位置出现与梯度成比例的不同长度(例如，下面的第一幅图)?

---

背景:

 #ggquiver示例图书馆(tidyverse)图书馆(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单位)来计算和绘制矢量字段.

---

结论:

要进行审查，我想获取任意的 rasterLayer 高程，将其转换为 data.frame ，计算 x 高度矢量字段的， y ， xmax 和 ymax 分量，这些分量调整箭头的大小，以便它们显示该点的相对斜率(例如在上面的图1和图2中)，并使用 ggquiver 进行绘制.像这样:

  names(r)<-"z"rd<-as.data.frame(r，xy = TRUE)#计算x，y，xend，yend为梯度向量，加到rd，然后绘制ggplot(rd)+geom_raster(aes(x，y，fill = z))+geom_quiver(aes(x，y，xend，yend)) 

解决方案

实际上，您要的是转换2D

现在要说明它可以在任意栅格上工作(假设已分配了投影)，让我们在此图像上对其进行测试，然后转换为栅格.这次，我们的分辨率较低，因此我们选择了较小的合计值.我们还将为最低值选择一种透明颜色，以绘制出更好的图:

  rast<-raster :: raster("https://i.stack.imgur.com/tXUXO.png")#添加一个假的任意投影，否则"terrain()"将不起作用:投影(rast)<--"+ proj = lcc + lat_1 = 48 + lat_2 = 33 + lon_0 = -100 + ellps = WGS84"raster2quiver(rast，聚合= 20，颜色= c(#FFFFFF00"，红色")) 

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

Problem statement:

With ggquiver::geom_quiver() we can plot vector fields, provided we know x, y, xend, and yend.

1. How can I calculate these parameters for an arbitrary RasterLayer of elevations?
2. How can I ensure that the size of these arrows indicates the slope for that particular vector such that the arrows appear different lengths proportional to the gradient at that location (e.g., the first plot below)?

---

Background:

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

A related appraoch uses rasterVis::vectorplot, which relies on raster::terrain (provided the field units == CRS units) to calculate and plot a vector field. Source code is here.

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

---

Conclusion:

To review, I'd like to take an arbitrary rasterLayer of elevation, convert it to a data.frame, calculate the x, y, xmax, and ymax components of an elevation vector field that size the arrows such that they show the relative slope at the point (as in plots 1 and 2 above), and plot with ggquiver. Something like:

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

解决方案

Effectively what you're asking is to convert a 2D scalar field into a vector field. There are a few different ways to do this.

The raster package contains the function terrain, which creates new raster layers that will give you both the angle of your desired vector at each point (i.e. the aspect), and its magnitude (the slope). We can use a little trigonometry to convert these into the North-South and East-West basis vectors used by ggquiver and add them to our original raster before turning the whole thing into a data frame.*

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)

However, in most cases this will not make for a good plot. If you don't first aggregate the raster, you will have one gradient for each pixel on the image, which won't plot well. On the other hand, if you do aggregate, you will have a nice vector field, but your raster will look "blocky". Therefore, having a single data frame for your plot is probably not the best way to go.

The following function will take a raster and plot it with an overlaid vector field. You can adjust how much the vector field is aggregated without affecting the raster, and you can specify an arbitrary vector of colours for your raster.

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

So, trying it out on your France example, after first defining a similar colour palette, we get

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

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

Now to show that it works on an arbitrary raster (provided it has a projection assigned) let's test it out on this image, as converted to a raster. This time, we have a lower resolution so we choose a smaller aggregate value. We'll also choose a transparent colour for the lowest values to give a nicer plot:

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

* _{I should point out that geom_quiver's mapping aesthetic takes arguments called u and v, which represent the basis vectors pointing North and East. The ggquiver package converts them to xend and yend values using stat_quiver. If you prefer to use xend and yend values you could just use geom_segment to plot your vector field, but this makes it more complex to control the appearance of the arrows. Hence, this solution will find the magnitude of the u and v values instead.}

这篇关于计算和绘制任意栅格图层的矢量场的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！