вычислить и построить векторное поле произвольного растрового слоя - PullRequest
Новая
       73

вычислить и построить векторное поле произвольного растрового слоя

7 голосов
/ 07 мая 2020

Формулировка проблемы:

С помощью 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()

enter image description here

A related appraoch uses rasterVis::vectorplot, which relies on raster::terrain (provided the field units == CRS units) to calculate and plot a vector field. Исходный код здесь . 

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

1 Ответ

5 голосов
/ 12 мая 2020

По сути, вы просите преобразовать 2D скалярное поле в векторное поле . Есть несколько разных способов сделать это.

Пакет растров содержит функцию 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)

Однако в большинстве случаев это не дает хорошего сюжета. Если вы сначала не объедините растр, у вас будет один градиент для каждого пикселя изображения, что не будет хорошо отображено. С другой стороны, если вы выполните агрегирование , у вас будет хорошее векторное поле, но ваш растр будет выглядеть «блочно». Следовательно, наличие одного фрейма данных для вашего графика, вероятно, не лучший способ go.

Следующая функция возьмет растр и построит его с наложенным векторным полем. Вы можете настроить степень агрегирования векторного поля, не влияя на растр, и вы можете указать произвольный вектор цветов для вашего растра. 

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 aestheti c принимает аргументы, называемые u и v, которые представляют базисные векторы, указывающие на север и восток. Пакет ggquiver преобразует их в значения xend и yend с использованием stat_quiver. Если вы предпочитаете использовать значения xend и yend, вы можете просто использовать geom_segment для построения векторного поля, но это усложняет управление внешним видом стрелок. Следовательно, это решение найдет вместо этого величину значений u и v.

...