Землекористування/покрив України у 2022 році

R
map
rayshader
visualisation
Візуалізація земляного покриву та землекористування в Україні
Автор

Ігор Мірошниченко

Дата публікації

2024-02-02

Важливо
  • Повний код для цієї візуалізації можна знайти на GitHub
  • Мапу у високій якості можна завантажити тут (51.8 Mb)
  • Натхнення брав у Milos Popovic

Карти землекористування та земляного покриву є важливим інструментом для вивчення взаємодії людини та природи. Вони дозволяють візуалізувати, як змінюється ландшафт внаслідок діяльності людини, а також виявляти тенденції у зміні землекористування. У цій статті я покажу, як можна побудувати вражаючі візуалізації земляного покриву та землекористування за допомогою R та пакету rayshader.

Пакети

Для цієї візуалізації я використовував наступні пакети:

if (!require("pacman")) install.packages("pacman")

pacman::p_load(
1    terra,
2    giscoR,
3    sf,
4    tidyverse,
5    ggtern,
6    elevatr,
7    png,
8    rayshader,
9    magick,
10    furrr,
11    here,
12    future,
13    rvest,
14    RSelenium,
15    netstat,
16    janitor,
17    extrafont
)
1
terra - для роботи з векторними та растровими даними
2
giscoR - для завантаження геоданих Євростату
3
sf - для роботи з векторними геоданими
4
tidyverse - для роботи з даними
5
ggtern - для побудови тернарних графіків
6
elevatr - для завантаження цифрових моделей рельєфу
7
png - для роботи з зображеннями
8
rayshader - для побудови 3D-візуалізацій
9
magick - для роботи з зображеннями
10
furrr - для паралельних обчислень
11
here - для роботи з шляхами
12
future - для паралельних обчислень
13
rvest - для парсингу веб-сторінок
14
RSelenium - для автоматизації браузера
15
netstat - для роботи з мережевими даними
16
janitor - для роботи з даними
17
extrafont - для роботи з шрифтами

Завантаження даних

Для візуалізації земляного покриву та землекористування я використовував дані Sentinel-2 10m land use/land cover time series of the world, які підготовлені Impact Observatory, Microsoft та Esri з використанням моделі глибокого навчання штучного інтелекту Impact Observatory для класифікації земель.

Кордони України

Для початку завантажимо кордони України за допомогою пакету giscoR:

country_sf <- gisco_get_countries(
1    country = "UA",
2    resolution = "1"
)

plot(st_geometry(country_sf))
1
country - код країни згідно з ISO 3166-1 alpha-2
2
resolution - рівень градації кордонів (1 - країна, 2 - область тощо)
Рисунок 1: Кордони України

Збережемо кордони України у файл ua-borders.png:

png("ua-borders.png")

Земляний покрив

Завантажимо растровий шар земляного покриву для України. Для цього необхідно перейти на сайт https://livingatlas.arcgis.com/landcoverexplorer та натиснути кнопку завантаження у лівому нижньому кутку карти.

Після чого відкриється мапа земляного покриву, на якій можна вибрати область для завантаження. Виберемо всі області, які входять до складу України. Клікаємо на кожну область щоб вибрати плитку та рік:

Копіюємо посилання та зберігаємо їх змінну urls:

urls <- c(
    "https://lulctimeseries.blob.core.windows.net/lulctimeseriesv003/lc2022/34U_20220101-20230101.tif",
    "https://lulctimeseries.blob.core.windows.net/lulctimeseriesv003/lc2022/35U_20220101-20230101.tif",
    "https://lulctimeseries.blob.core.windows.net/lulctimeseriesv003/lc2022/36U_20220101-20230101.tif",
    "https://lulctimeseries.blob.core.windows.net/lulctimeseriesv003/lc2022/37U_20220101-20230101.tif",
    "https://lulctimeseries.blob.core.windows.net/lulctimeseriesv003/lc2022/35T_20220101-20230101.tif",
    "https://lulctimeseries.blob.core.windows.net/lulctimeseriesv003/lc2022/36T_20220101-20230101.tif",
    "https://lulctimeseries.blob.core.windows.net/lulctimeseriesv003/lc2022/37T_20220101-20230101.tif"
)

Для завантаження даних використаємо функцію download.file(). За замовчування вона завантажує файли впродовж 60 секунд і якщо файл не встиг завантажитись, видасть помилку. Щоб уникнути цього, збільшимо час завантаження до 480 секунд:

options(timeout = 480)

Створимо функцію для завантаження растрових даних якщо вони ще не були завантажені. Це дозволить нам уникнути повторного завантаження файлів, якщо вони вже були завантажені. Імена файлів будуть взяті з посилань за допомогою функції basename():

download_if_not_exists <- function(url) {
    if (!file.exists(basename(url))) {
        download.file(url, basename(url), mode = "wb")
    }
}

Завантажимо файли за допомогою нашої функції та функції map() з пакету purrr:

map(urls, download_if_not_exists)

Об’єднання растрових файлів

Збережемо перелік всіх завантажених файлів у змінну raster_files:

raster_files <- list.files(
1    path = here(),
2    pattern = "20230101.tif$",
3    full.names = TRUE
)
1
path - шлях до папки, в якій знаходяться файли. В даному випадку це поточна папка, яку ми отримуємо за допомогою функції here()
2
pattern - шаблон імен файлів, які ми шукаємо. В даному випадку це файли, які закінчуються на 20230101.tif
3
full.names - будемо зберігати повні шляхи до файлів

Для подальшої роботи і об’єднання всіх файлів в один створимо змінну з crs в якій буде збережено інформацію про систему координат, в якій знаходяться файли. В даному випадку це EPSG:32635:

crs <- "EPSG:4326"

Об’єднаємо всі файли в один за допомогою функції terra::rast():

if(!file.exists("ua_land_cover_vrt.vrt")){
    for(raster in raster_files){
        cat("Start with:", raster, "\n")
1        rasters <- rast(raster)
2        country <- country_sf %>%
            sf::st_transform(
                crs = crs(
                    rasters
                )
            )
3        land_cover <- crop(
            rasters,
            vect(
                country
            ),
            snap = "in",
            mask = TRUE
        ) %>% 
4        terra::aggregate(
            fact = 5,
            fun = "modal"
        ) %>% 
5        terra::project(crs)
        # write raster
6        terra::writeRaster(
            land_cover,
            paste0(
                raster,
                "_ua",
                ".tif"
            ),
            overwrite=FALSE
        )
    }
  }
1
Завантажуємо растровий файл за допомогою функції terra::rast()
2
Трансформуємо систему координат.
3
Обрізаємо растровий файл за допомогою векторного файлу кордонів України та зберігаємо його.
4
Агрегуємо растровий файл для зменшення об’єму даних.
5
Проектуємо растровий файл в систему координат EPSG:4326.
6
Зберігаємо растровий файл.

Нарешті можемо об’єднати всі файли в один за допомогою функції terra::vrt(), яка створить віртуальний растровий файл:

r_list <- list.files(
    path = here(),
    pattern = "_ua.tif",
    full.names = TRUE
)

land_cover_vrt <- terra::vrt(
    r_list,
    "ua_land_cover_vrt.vrt",
    overwrite = TRUE
)

Отримання оригінальної палітри кольорів

Для візуалізації земляного покриву нам знадобиться палітра кольорів, яка використовується на оригінальній мапі. Для цього завантажимо один з растрових файлів та витягнемо з нього палітру кольорів:

1ras <- rast(
    raster_files[[1]]
)

2raster_color_table <- do.call(
    data.frame,
3    coltab(ras)
)
1
Завантажуємо растровий файл за допомогою функції terra::rast()
2
Створюємо датафрейм з палітрою кольорів
3
Витягуємо палітру кольорів з растрового файлу

Зверніть увагу на те, що палітра кольорів зберігається у форматі RGB:

head(raster_color_table)

Для подальшої роботи нам знадобиться HEX-код кольорів, тому перетворимо RGB-коди кольорів у HEX-коди за допомогою функції ggtern::rgb2hex():

hex_code <- rgb2hex(
    r = raster_color_table[,2],
    g = raster_color_table[,3],
    b = raster_color_table[,4]
)

hex_code[1:14]

Як бачимо, тільки перші 11 спостережень мають HEX-коди. Отже нам потрібні тільки 2-3, 4-5 та 7-11 рядки. Збережемо їх у змінну cols:

cols <- hex_code[c(2:3, 5:6, 8:12)]

Тепер ми можемо використати ці кольори для візуалізації земляного покриву:

1from <- c(1:2, 4:5, 7:11)
2to <- t(col2rgb(cols))

3land_cover_vrt <- na.omit(land_cover_vrt)

4land_cover_ua <- terra::subst(
    land_cover_vrt,
    from = from,
    to = to,
    names = cols
)
1
Вибираємо рядки з HEX-кодами кольорів
2
Перетворюємо HEX-коди кольорів у RGB-коди
3
Видаляємо пропущені значення
4
Замінюємо кольори у растровому файлі

Тепер можемо поглянути на отриманий растровий файл:

terra::plotRGB(land_cover_ua)

Цифрова модель рельєфу

Для побудови 3D-візуалізації нам знадобиться цифрова модель рельєфу. Завантажимо її за допомогою функції elevatr::get_elev_raster():

elev <- elevatr::get_elev_raster(
1    locations = country_sf,
2    z = 8,
3    clip = "locations"
)
1
locations - векторний файл з кордонами України
2
z - масштаб цифрової моделі рельєфу. Для великих територій рекомендується зменшувати масштаб.
3
clip - обрізка за кордонами України

Для об’єднання растрових файлів земляного покриву та цифрової моделі рельєфу нам необхідно перетворити у відповідності один до одного. Для цього ми використаємо проекцію, яку можна знайти за посиланням https://epsg.io/ для своєї країни:

crs_lambert <-
    "+proj=tmerc +lat_0=0 +lon_0=21 +k=1 +x_0=300000 +y_0=0 +ellps=krass +towgs84=24,-121,-76,0,0,0,0 +units=m +no_defs +type=crs"

Тепер ми можемо перетворити цифрову модель рельєфу у відповідність до растрового файлу земляного покриву:

land_cover_ua_resampled <- terra::resample(
1    x = land_cover_ua,
2    y = terra::rast(elev),
3    method = "near"
) %>% 
4    terra::project(crs_lambert)
1
Растровий файл земляного покриву
2
Цифрова модель рельєфу
3
Метод перетворення
4
Проекція

Подивимось на результат:

terra::plotRGB(land_cover_ua_resampled)

Збережемо рисунок у файл land_cover_ua.png і зчитаємо його за допомогою пакету png:

img_file <- "land_cover_ua.png"

terra::writeRaster(
    land_cover_ua_resampled,
    img_file,
    overwrite = TRUE,
    NAflag = 255
)

img <- png::readPNG(img_file)

Створення 3D-візуалізації

Перш ніж створити 3D-візуалізацію, нам необхідно переконатися, що наша растрова модель має однакову проекцію з цифровою моделлю рельєфу. Також нам потрібно перетворити растрову модель у матрицю, щоб пакет rayshader міг працювати з нею:

elev_lambert <- elev %>%
   terra::rast() %>%
   terra::project(crs_lambert)

elmat <- rayshader::raster_to_matrix(
   elev_lambert 
)

Визначимо висоту та ширину растрової моделі:

h <- nrow(elev_lambert)
w <- ncol(elev_lambert)

Тепер ми можемо створити 3D-візуалізацію за допомогою пакету rayshader:

elmat %>% 
    rayshader::height_shade(
        texture = colorRampPalette(
            cols[9]
        )(256)
    ) %>% 
    rayshader::add_overlay(
        img,
        alphalayer = 1
    ) %>% 
    rayshader::plot_3d(
        elmat,
        zscale = 12,
        solid = FALSE,
        shadow = TRUE,
        shadow_darkness = 1,
        background = "white",
        windowsize = c(
            w / 10, h / 10
        ),
        zoom = .5,
        phi = 85,
        theta = 0
    )

rayshader::render_camera(
    zoom = .58
)

В результаті ми отримаємо 3D-візуалізацію земляного покриву та цифрової моделі рельєфу.

Рендер рисунку

Тепер ми готові зберегти все у високоякісний рисунок. Збережемо його у файл 3d_land_cover_ua-dark.png за допомогою функції rayshader::render_highquality() і використаємо для цього високоякісну текстуру освітлення, яка збережена у файлі air_museum_playground_4k.hdr:

filename <- "3d_land_cover_ua-dark.png"

rayshader::render_highquality(
    filename = filename,
    preview = TRUE,
    light = FALSE,
    environment_light = here('air_museum_playground_4k.hdr'),
    intensity_env = 1,
    rotate_env = 90,
    interactive = FALSE,
    parallel = TRUE,
    width = w * 1.5,
    height = h * 1.5
)

Збір статистичної інформації

Парсинг

Також хотілося б зібрати статистичну інформацію про земляний покрив України щоб у подальшому додати її у легенду.

Для цього зберемо інформацію зі сторінки опису даних на сайті Sentinel-2 10m Land Use/Land Cover Time Series.

Сайт динамічно підтягує контент, тому нам знадобиться пакет RSelenium для парсингу даних. Відкриємо вікно браузера та перейдемо на сайт:

driver <- rsDriver(browser = "firefox",
                   chromever = "114.0.5735.90",
                   verbose = FALSE,
                   port = free_port())

remDr <- driver$client

remDr$navigate("https://www.arcgis.com/home/item.html?id=cfcb7609de5f478eb7666240902d4d3d")

Тепер можемо знайти елементи сторінки за допомогою CSS-селекторів та зібрати інформацію з таблиці:

data_table <- remDr$findElement(using = "xpath", '/html/body/div[3]/div/div[2]/div/div[2]/div/main/div[2]/div[2]/div[1]/div/div/div/div[8]/table/tbody')

Тепер можемо зібрати інформацію з таблиці:

class_tbl <- data_table$getPageSource() %>% 
  unlist() %>%
  read_html() %>% 
  html_table() %>% 
  .[[1]] %>%
  row_to_names(row_number = 1) %>% 
  mutate(Value = as.factor(Value))

class_tbl
# A tibble: 9 × 3
  Value Name               Description                                          
  <dbl> <chr>              <chr>                                                
1     1 Water              Areas where water was predominantly present througho…
2     2 Trees              Any significant clustering of tall (~15 feet or high…
3     4 Flooded vegetation Areas of any type of vegetation with obvious intermi…
4     5 Crops              Human planted/plotted cereals, grasses, and crops no…
5     7 Built Area         Human made structures; major road and rail networks;…
6     8 Bare ground        Areas of rock or soil with very sparse to no vegetat…
7     9 Snow/Ice           Large homogenous areas of permanent snow or ice, typ…
8    10 Clouds             No land cover information due to persistent cloud co…
9    11 Rangeland          Open areas covered in homogenous grasses with little…

Зверніть увагу, що змінна Value не містить записів 3 та 6.

Підрахунок площі

Тепер можемо підрахувати площу кожного класу земляного покриву. На випадок, якщо ми захочемо порахувати площі окремо для кожної області, додамо id до змінної country_sf:

country_sf$id <- 1:nrow(country_sf)

Зберемо всі класи земляного покриву у змінну classes:

classes <- class_tbl$Value

Проведемо підрахунок статистики для кожного класу земляного покриву за допомогою функції exactextractr::exact_extract():

zonal_stats_ukr <- exact_extract(
  land_cover_vrt,
  country_sf
)

Тепер можемо підрахувати пропорції площі кожного класу земляного покриву. Для цього використаємо паралельні обчислення за допомогою пакету future:

plan(multisession)

ukraine_multiclass <- future_map_dfr(zonal_stats_ukr, function(x) {
  as.data.frame(
    prop.table(
      table(factor(x[, 1], levels = classes))
    )
  )
})

Об’єднаємо результати з попередньою таблицею, додамо змінну colors, яка буде містити палітру кольорів для кожного класу земляного покриву, а також змінну perc, яка буде містити відсотове співвідношення площі кожного класу земляного покриву. Відсортуємо результати за спаданням площі та видалимо класи, які не мають значення:

ukraine_multiclass_sf <- ukraine_multiclass %>% 
  as_tibble() %>%
  left_join(
    class_tbl,
    by = c("Var1" = "Value")
  ) %>% 
  mutate(colors = c(
    "#419bdf", "#397d49", "#7a87c6", 
    "#e49635", "#c4281b", "#a59b8f", 
    "#a8ebff", "#616161", "#e3e2c3"
  ), .after = Name,
  perc = scales::percent(Freq / sum(Freq), accuracy = .01, trim = FALSE)) %>% 
  arrange(desc(Freq)) %>% 
  filter(!(Name %in% c("Bare ground", 'Snow/Ice', 'Clouds')))

ukraine_multiclass_sf
# A tibble: 6 × 6
   Var1    Freq Name               colors  perc   Description                   
  <dbl>   <dbl> <chr>              <chr>   <chr>  <chr>                         
1     5 0.665   Crops              #e49635 66.45% Human planted/plotted cereals…
2     2 0.189   Trees              #397d49 18.91% Any significant clustering of…
3    11 0.0706  Rangeland          #e3e2c3 7.06%  Open areas covered in homogen…
4     7 0.0504  Built Area         #c4281b 5.04%  Human made structures; major …
5     1 0.0212  Water              #419bdf 2.12%  Areas where water was predomi…
6     4 0.00392 Flooded vegetation #7a87c6 0.39%  Areas of any type of vegetati…

Збираємо фінальну візуалізацію

Легенда

Тепер можемо побудувати легенду для нашої візуалізації. Для цього використаємо пакет ggplot2:

plot <- ukraine_multiclass_sf %>% 
  ggplot(aes(x = Freq, y = fct_reorder(Name, Freq))) +
  geom_col(aes(fill = fct_reorder(Name, Freq))) +
  geom_text(
    aes(label = perc), size = 8, hjust = -0.1, vjust = 0.5, family = "Fira Sans", fontface = "bold") +
  scale_fill_manual(values = rev(ukraine_multiclass_sf$colors)) +
  theme_void() +
  coord_cartesian(clip = "off") +
  scale_x_continuous(expand = c(.01, .01)) +
  theme(
    legend.position = "none",
    axis.text.y = element_text(size = 28, hjust = 1, family = "Georgia"),
    plot.margin = margin(15, 85, 15, 15))

plot

Збережемо легенду у файл land_cover_legend.png:

legend_name <- "land_cover_legend.png"
ggsave(legend_name, plot, width = 10, height = 6)

Збірка фінальної візуалізації

За допомогою пакету magick можемо зберегти всі елементи візуалізації у один файл:

lc_img <- magick::image_read(
  filename
)

my_legend <- magick::image_read(
  legend_name
)

my_legend_scaled <- magick::image_scale(
  magick::image_background(
    my_legend, "none"
  ), 2500
)

p <- magick::image_composite(
    magick::image_scale(
        lc_img, "x7000" 
    ),
    my_legend_scaled,
    gravity = "southwest",
    offset = "+1000+1000"
    ) %>% 
    image_annotate(
    "Land cover in 2022",
    size = 150,
    color = alpha("#e49635", .5),
    font = "Georgia",
    gravity = "north",
    location = "+0+300"
  ) %>% 
  image_annotate(
    "Ukraine",
    size = 300,
    color = "#e49635",
    font = "Georgia",
    gravity = "north",
    location = "+0+500"
  ) %>% 
  image_annotate(
    "©2023 Ihor Miroshnychenko (https://aranaur.rbind.io)",
    size = 100,
    color = alpha("grey20", .75),
    font = "Georgia",
    gravity = "southeast",
    location = "+100+200"
  ) %>% 
  image_annotate(
    "Data: Esri | Sentinel-2 Land Cover Explorer",
    size = 100,
    color = alpha("grey20", .75),
    font = "Georgia",
    gravity = "southeast",
    location = "+100+100"
  )

Збережемо фінальну візуалізацію у файл 3d_ua_land_cover_final_eng.png:

magick::image_write(
  p, "3d_ua_land_cover_final_eng.png"
)

В результаті ми отримаємо фінальну візуалізацію земляного покриву України:

Цитата

BibTeX citation:
@misc{miroshnychenko2024,
  author = {Miroshnychenko, Ihor},
  title = {Ukraine’s land use/cover in 2022},
  date = {2024-02-02},
  url = {https://aranaur.rbind.io/blog/2024/01},
  langid = {uk-UA}
}
Будь-ласка, цитуйте цю роботу як:
Miroshnychenko, Ihor. 2024. “Ukraine’s land use/cover in 2022.” February 2, 2024. https://aranaur.rbind.io/blog/2024/01.