Визуализация трехмерных моделей земной поверхности

Примеры 3D-визуализаций цифровой модели рельефа местности средствами библиотеки rayvista

rayvista
spatial
Rstats
Автор

Е.Н. Матеров

Дата публикации

4 марта 2024 г.

Библиотека rayvista

В R есть множество совершенно изумительных библиотек с открытым кодом, одна из них – rayshader, которая позволяет делать потрясающие 3D-изображения на основе матрицы высот и различных алгоритмов отмывки теней и трассировки лучей. Помимо создания 3D-картографических моделей, в rayshader можно преобразовывать объекты ggplot2 в красивые трехмерные визуализации данных. Полученные модели можно вращать и исследовать в интерактивном режиме, а движение камеры можно запрограммировать для создания анимации. Подробно работу с rayshader, функции, руководства и результаты работы можно посмотреть на сайте библиотеки, на персональном сайте автора библиотеки Rayverse Blog и в главе Трехмерные модели интерактивной книги Визуализация и анализ географических данных на языке R.

Библиотека rayvista создана как плагин к rayshader и содержит основную функцию plot_3d_vista() для создания 3D-визуализации практически любого местоположения на Земле по географическим координатам. Так были созданы альбомы Горы России и Непокоренные вершины мира. Для создания цифровой модели рельефа местности на основе высот над уровнем моря в rayvista используется библиотека elevatr. Для создания карт, которые накладываются на 3D-основу высот, используется библиотека maptiles, которая позволяет загружать, комбинировать и отображать карты различных провайдеров (по умолчанию это Esri, а также OpenStreetMap, CARTO и Thunderforest).

Полученные модели можно комбинировать с командами rayshader, накладывать объекты OpenStreetMap, и т.д. Установить библиотеку можно с помощью команды

# install.packages("devtools")
devtools::install_github("h-a-graham/rayvista", dependencies=TRUE)

При создании модели открывается RGL-окно, например, в macOS необходимо установить XQuartz.

Базовая работа с 3D-моделями

Создадим простейшую 3D-модель. В качестве исходной визуализации выберем модель вершины Амадаблам – вершины в Гималаях, высота главного пика которой равна 6 814 м.

Рисунок 1: Пример 3D-модели в rayvista – вершина Амадаблам 
1Amadablam <- plot_3d_vista(lat = 27.8599,
                           long = 86.8614, 
2                           zscale = 6,
3                           zoom = 0.6,
                           soliddepth = 4000,
4                           radius = 5000,
5                           overlay_detail = 14,
6                           theta = -65, phi = 25,
7                           windowsize = 1200,
                           background = "grey10")
1
широта и долгота, выраженная в WGS84;
2
при условии, что x = y, zscale показывает соотношение между x и z (либо y и z), по умолчанию zscale = 2;
3
коэффициент масштабирования, по умолчанию равен 1;
4
радиус, определяющий ограничивающую область;
5
число между 0 и 20, отвечающее за масштаб в maptiles::get_tiles(), уровни масштабирования можно посмотреть в OpenStreetMap wiki, по умолчанию значение равно 13;
6
угол вращения и азимутальный угол соответственно;
7
положение, ширина и высота устройства RGL, отображающего модель.

В библиотеке rayvista возможно использовать карты-подложки различных провайдеров (например, OpenStreetMap, Esri, CARTO или Thunderforest, также, см. статью по поводу доступности карт Stamen). Все возможные слои показывает команда

?maptiles::get_tiles

В примере ниже покажем различные варианты для визуализации острова Татышев, г. Красноярск.

Krasnoyarsk <- plot_3d_vista(lat = 56.01986, 
                             long = 92.932413, 
                             zscale = 3, 
                             zoom = 0.5,
1                             img_provider = "Thunderforest.Landscape",
2                             api_key = "your_api_key_value",
                             overlay_detail = 14,
                             theta = -40,  
                             phi = 25,
                             windowsize = 1300,
                             soliddepth = 50,
3                             # outlier_filter = 0.0001,
4                             fill_holes = TRUE)
1
вид картографической подолжки;
2
значение API-ключа для Thunderforest-карт, можно получить на сайте, для других карт API-ключ не требуется;
3
если в модели появляются артефакты, это значит что возвращаемые данные о рельефе имеют ошибочно низкие значения, тогда нужно выставить значение от 0 до 1, ближе к 0;
4
заполняет значения NA полученные в цифровой модели рельефа местности.
Рисунок 2: Пример 3D-модели с картой Thunderforest.Landscape, где показаны топографические особенности

Отметим, что если в предыдущем коде, который формирует 3D-изображение, добавить show_vista = FALSE (в команде plot_3d_vista()), то мы не увидим RGL-окно, а сохраним объект, у которого будут атрибуты dem_matrix и texture. Теперь восстановить изображение на рисунке (см. рис. 2) можно так:

Krasnoyarsk$dem_matrix |>
  height_shade() %>% 
  add_overlay(., Krasnoyarsk$texture,
              rescale_original = TRUE) %>%
  plot_3d(., Krasnoyarsk$dem_matrix, zscale = 3,
          windowsize = 1300, 
          soliddepth = 50,
          zoom = 0.5, phi = 25, theta = -40)

Сравним предыдущую карту с другими картами-подложками.

Рисунок 3: Пример 3D-модели с картой Thunderforest.Outdoors
Рисунок 4: Пример 3D-модели с картой OpenTopoMap
Рисунок 5: Пример 3D-модели с картой OpenStreetMap

Для лучшего восприятия можно оставить только 3D-поверхность модели, указав solid = FALSE как в следующем примере.

Код 
MalySemyachik <- plot_3d_vista(lat  = 54.119, 
                               long = 159.656, 
                               zscale = 3, 
                               zoom = 0.6,
                               radius = 5000,
                               overlay_detail = 14,
                               theta = -50,
                               phi = 30,
                               windowsize = 1100, 
                               solid = FALSE,
                               fill_holes = TRUE)
Рисунок 6: Стратовулкан Малый Семячик с горным озером на Камчатке

Для имитации эффекта размытия можно использовать функцию render_depth().

Код 
HopkinsNZ <- plot_3d_vista(lat = -44.042238, 
                           long = 169.860985, 
                           radius = 5000, 
                           overlay_detail = 14,
                           elevation_detail = 13, 
                           zscale = 5, 
                           theta = 25, phi=25, zoom = 0.6,
                           windowsize = 1200, 
                           solid = T, 
                           background = 'white')

render_depth(focus = 0.6, focallength = 15, clear=TRUE)
Рисунок 7: Mount Hopkins (Новая Зеландия)

Нанесение дополнительных географических данных на 3D-визуализации

Поднятие плоскости воды

Одна из интересных возможностей, для которой может быть полезна библиотека rayvista – это оценка природных рисков наводнений локальных территорий путем параллельного поднятия плоскости воды. В качестве примера рассмотрим г. Улан-Удэ. Выберем точку внутри города с определенными географическими координатами и найдем оценочную высоту над уровнем моря с помощью сервиса Elevation Finder, в нашем случае это 490 м (в rayshader также существует команда detect_water()). Получим модель местности и поднимем уровень на 20 м до 510 м, чтобы сделать оценку риска затопления.

Рисунок 8: Оценка риска затопления с помощью 3D-моделирования для г. Улан-Удэ увеличением уровня воды на +20 м, возможный уровень затопления показан красным цветом 
UlanUde <- plot_3d_vista(lat  = 51.834, 
                         long = 107.5696, 
                         zscale = 2, 
                         zoom = 0.5,
                         overlay_detail = 14,
                         theta = -65, 
                         radius = 5500,
                         img_provider = "Thunderforest.Outdoors",
                         api_key = "your_api_key",
                         windowsize = 1400, 
                         soliddepth = 300,
                         phi = 25,
1                         water = TRUE,
2                         waterdepth = 510,
3                         wateralpha = 0.3,
4                         watercolor = "red",
5                         waterlinecolor = "white",
                         fill_holes = TRUE)
1
добавление слоя воды;
2
уровень воды;
3
прозрачность слоя воды;
4
цвет слоя воды;
5
цвет линий по краям слоя воды.

Код ниже добавляет маркеры и метки к текущему 3D-образу.

Код 
render_label(heightmap = UlanUde,
             lat  = 51.824, 
             long = 107.525,
             altitude = 140, 
             text = "Улан-Удэ (+20 метров)",
             extent = attr(UlanUde, 'extent'),
             textsize = 1.6, linewidth = 4,
             textcolor = "grey20", linecolor = "grey20")

Нанеcение данных OpenStreetMap

На 3D-образы можно наносить помимо маркеров, например, растровые данные или данные из OpenStreetMap.

Предостережение

Практически все сервисы стандарта Web Map Tile Service, такие как OSM, Stamen, GoogleMaps и т.д., предоставляют тайлы в системе координат (CRS) EPSG:3857 (WebMercator). Это означает, что если объект нового слоя находится в другой проекции CRS, который не является EPSG:3857 (например, EPSG:4326), необходимо привести наносимый объект (sf или растр) и подложку к единой проекции путем перепроектирования, иначе модель будет будет выглядеть деформированной.

В качестве примера покажем, как сделать дополнительный слой на 3D-модели, который представляет собой данные по дорогам, загруженным из сервиса OpenStreetMap. Выберем Сан-Франциско (Калифорния, США), для которого характерен сложный рельеф местности; по некоторым сведениям в Сан-Франциско насчитывают около пятидесяти возвышенностей.

Загрузим 3D-модель.

Рисунок 9: Модель Сан-Франциско 
SanFrancisco <- plot_3d_vista(lat = 37.75, 
                              long = -122.45,
                              zscale = 3, 
                              zoom = 0.5,
                              img_provider = "Thunderforest.Outdoors",
                              api_key = "your_api_key",
                              overlay_detail = 14,
                              theta = -40, 
                              phi = 25,
                              windowsize = 1300,
                              show_vista = FALSE)

Определим прямоугольник, ограничивающий область.

library(sf)

bbox <- st_bbox(c(xmin = -122.5129, 
                  ymin = 37.70026, 
                  xmax = -122.3871, 
                  ymax = 37.79975),  
                crs = st_crs(4326))

Теперь загрузим дорожную сеть из OpenStreetMap.

library(osmdata)

SanFrancisco_highway <- opq(bbox) |>
  add_osm_feature("highway") |>
  osmdata_sf()

Следующий этап – определение новых проекций, без этого график будет искажен.

SanFrancisco_lines = 
  st_transform(SanFrancisco_highway$osm_lines, 3857)

bbox_transformed = st_bbox(
  st_transform(
    st_as_sfc(bbox), 
    3857
  )
)

librarary(raster)

extent_zoomed = extent(bbox_transformed)

Окончательно:

SanFrancisco$dem_matrix |>
  height_shade() %>% 
  add_overlay(., SanFrancisco$texture,
              rescale_original = TRUE) %>%
  # дороги
  add_overlay(
    generate_line_overlay(
      SanFrancisco_lines, 
      extent = extent_zoomed,
      linewidth = 1.2, color = "black",
      heightmap = SanFrancisco$dem_matrix)) %>%
  plot_3d(., SanFrancisco$dem_matrix, 
          zscale = 3,
          windowsize = 1300, 
          zoom = 0.5, phi = 25, theta = -40)

Теперь можно сравнить модель после нанесения дорожной сети на образ.

Рисунок 10: Модель Сан-Франциско с нанесенными OSM-дорогами

Нанесение дополнительных слоев на 3D-модели

Несомненно, на модели можно накладывать и другие слои с исходными данными или результатами моделирования, что делает технологию весьма привлекательной (см., например, статью Adding Open Street Map Data to Rayshader Maps). Приведем еще один пример, сделаем слой на 3D-модели, показывающий происшествия, на которые отреагировала пожарная служба Сан-Франциско. Отметим, соответствующий набор данных также содержит географические координаты.

data_sf |>
  dplyr::select(`Incident Date`, point)
# A tibble: 616,630 × 2
   `Incident Date`     point                            
   <dttm>              <chr>                            
 1 2008-04-01 00:00:00 POINT (-122.41837339 37.74208979)
 2 2008-04-01 00:00:00 POINT (-122.39489 37.756291)     
 3 2008-04-01 00:00:00 POINT (-122.407468 37.78008)     
 4 2008-04-01 00:00:00 POINT (-122.42684908 37.77612642)
 5 2008-04-01 00:00:00 POINT (-122.4863941 37.77428492) 
 6 2008-04-01 00:00:00 POINT (-122.4481912 37.7597267)  
 7 2008-04-01 00:00:00 POINT (-122.405223 37.788694)    
 8 2008-04-01 00:00:00 POINT (-122.416457 37.739056)    
 9 2008-04-01 00:00:00 POINT (-122.392082 37.781846)    
10 2008-04-01 00:00:00 POINT (-122.46678342 37.75300898)
# ℹ 616,620 more rows

Сделаем небольшую предобработку данных.

library(magrittr)

data_sf$coordinates <- 
  str_extract(data_sf$point, 
    "(?<=\\()([^()]*?)(?=\\)[^()]*$)")

data_sf <- data_sf |> 
  separate_wider_delim(coordinates, " ", 
    names = c("lon", "lat"))

data_sf$lat %<>% as.numeric()
data_sf$lon %<>% as.numeric()

В качестве примера, выделим точки вызовов, отвечающие 2022 году.

points_sf <- data_sf |>
  dplyr::select(`Incident Date`, lat, lon) |>
  mutate(year = lubridate::year(`Incident Date`)) |> 
  dplyr::filter(year == 2022) |>
  na.omit() |>
  st_as_sf(coords = c("lon", "lat"), crs = 4326) %>% 
  st_transform(., 3857)
points_sf
Simple feature collection with 33260 features and 2 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -13685560 ymin: 0 xmax: 13627530 ymax: 4555811
Projected CRS: WGS 84 / Pseudo-Mercator
# A tibble: 33,260 × 3
   `Incident Date`      year            geometry
 * <dttm>              <dbl>         <POINT [m]>
 1 2022-01-16 00:00:00  2022 (-13625484 4548073)
 2 2022-01-16 00:00:00  2022 (-13626605 4547332)
 3 2022-01-21 00:00:00  2022 (-13625214 4547874)
 4 2022-01-21 00:00:00  2022 (-13629079 4549058)
 5 2022-01-23 00:00:00  2022 (-13628711 4551584)
 6 2022-01-30 00:00:00  2022 (-13627494 4538682)
 7 2022-02-08 00:00:00  2022 (-13625338 4547751)
 8 2022-02-08 00:00:00  2022 (-13627882 4549346)
 9 2022-02-08 00:00:00  2022 (-13627827 4546331)
10 2022-02-10 00:00:00  2022 (-13627366 4548730)
# ℹ 33,250 more rows
SanFrancisco$dem_matrix |>
  height_shade() %>% 
  add_overlay(., SanFrancisco$texture,
              rescale_original = TRUE) %>%
  # пожары
  add_overlay(
    generate_point_overlay(
      points_sf, extent = extent_zoomed,
      size = 0.8, color = "#C11317",
      heightmap = SanFrancisco$dem_matrix)) %>%
  plot_3d(., SanFrancisco$dem_matrix, 
          zscale = 3,
          windowsize = 1300, 
          zoom = 0.5, phi = 25, theta = -40)
Рисунок 11: 3D-модель с нанесенными происшествиями, на которые отреагировала пожарная служба Сан-Франциско в 2022 году

Как видно на рисунке (см. рис. 11), наибольшее количество происшествий приходится на деловой центр северо-востока города.

Заключение

Библиотека rayvista позволяет довольно легко формировать 3D-визуализацию для различных мест планеты Земля и, при необходимости, моделировать различные физические процессы, что делает ее весьма перспективным решением. Мы рассмотрели некоторые примеры применения библиотеки, дополнительные возможности и нанесение дополнительной географической информации на слои.

Хорошим руководством по rayvista может также стать видео 3D digital elevation maps with elevatr, rayshader and rayvista in R, в котором Milos Popovic рассказывает этапы создания реалистичных цифровых 3D-карт высот любого местоположения, области или страны в мире с использованием библиотек elevatr, rayshader и rayvista в R.

─ Session info ───────────────────────────────────────────────────────────────
 setting  value
 version  R version 4.3.2 (2023-10-31)
 os       macOS Monterey 12.1
 system   aarch64, darwin20
 ui       X11
 language (EN)
 collate  ru_RU.UTF-8
 ctype    ru_RU.UTF-8
 tz       Asia/Krasnoyarsk
 date     2024-03-01
 pandoc   2.18 @ /Users/materov/opt/miniconda3/envs/ox/bin/ (via rmarkdown)
 quarto   1.4.550

─ Packages ───────────────────────────────────────────────────────────────────
 package     * version date (UTC) lib source
 dplyr       * 1.1.4   2023-11-17 [1] CRAN (R 4.3.1)
 forcats     * 1.0.0   2023-01-29 [1] CRAN (R 4.3.0)
 ggplot2     * 3.5.0   2024-02-23 [1] CRAN (R 4.3.1)
 lubridate   * 1.9.3   2023-09-27 [1] CRAN (R 4.3.1)
 magrittr    * 2.0.3   2022-03-30 [1] CRAN (R 4.3.0)
 purrr       * 1.0.2   2023-08-10 [1] CRAN (R 4.3.0)
 readr       * 2.1.5   2024-01-10 [1] CRAN (R 4.3.1)
 sessioninfo * 1.2.2   2021-12-06 [1] CRAN (R 4.3.0)
 sf          * 1.0-15  2023-12-18 [1] CRAN (R 4.3.1)
 stringr     * 1.5.1   2023-11-14 [1] CRAN (R 4.3.1)
 tibble      * 3.2.1   2023-03-20 [1] CRAN (R 4.3.0)
 tidyr       * 1.3.1   2024-01-24 [1] CRAN (R 4.3.1)
 tidyverse   * 2.0.0   2023-02-22 [1] CRAN (R 4.3.0)

 [1] /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library

──────────────────────────────────────────────────────────────────────────────

Ссылка для цитирования

BibTeX
@online{матеров2024,
  author = {Матеров, Е.Н.},
  title = {Визуализация трехмерных моделей земной поверхности},
  date = {2024-02-02},
  url = {https://www.naukaidannye.netlify.app/blog/posts/2024-02-02-rayvista},
  langid = {ru}
}
На публикацию можно сослаться как
Матеров Е.Н. Визуализация трехмерных моделей земной поверхности [Electronic resource]. 2024. URL: https://www.naukaidannye.netlify.app/blog/posts/2024-02-02-rayvista.