Код
# общие библиотеки
library(tidyverse)
library(magrittr)
# работа внутри директорий
library(here)
# библиотеки для работы с географическими данными
library(raster)
library(marmap)
library(elevatr)
library(sp)
library(sf)
library(terra)Исследование повторяемости ветров в южной части Российской Федерации
Замечание
Данный документ создан на языке программирования R с помощью воспроизводимых вычислений и опубликован на основе издательской системы Quarto версии 1.3.228 на сервисе Netlify. При изменении исходных данных, документ автоматически будет сгенерирован с учетом изменений.
Исходные данные
Загрузим необходимые библиотеки языка программирования R.
Загрузим табличные данные.
Код
data_path <- here("data")
data_table_path <- here(data_path, "wind_2021.csv")
wind_2021_raw <- read_csv(data_table_path)
wind_2021 <- wind_2021_raw
# правильные типы данных
wind_2021$lon %<>% as.numeric()
wind_2021$lat %<>% as.numeric()
# месяцы как факторы (упорядочивание)
wind_2021$month <- factor(wind_2021$month,
labels = c(
"апрель",
"май",
"июнь",
"июль"
),
levels = c(
"апрель",
"май",
"июнь",
"июль"
)
)Посмотрим на структуру исходных табличных данных. Здесь:
| переменная | значение |
|---|---|
lat |
широта точки |
lon |
долгота |
wind_value |
повторяемость ветра в 2021 г. (количество часов) |
year |
год |
month |
месяц |
Данные представляют собой объединенную таблицу, в которой каждая строка соответствует одному наблюдению, а названия столбцов – переменным.1 Такие данные удобно преобразовывать и использовать в дальнейшем для графиков.
Код
wind_2021# A tibble: 90,720 × 5
lat lon wind_value year month
<dbl> <dbl> <dbl> <dbl> <fct>
1 60 45 1 2021 апрель
2 60 45.2 52.4 2021 апрель
3 60 45.5 102. 2021 апрель
4 60 45.8 153 2021 апрель
5 60 46 194. 2021 апрель
6 60 46.2 197. 2021 апрель
7 60 46.5 155. 2021 апрель
8 60 46.8 124. 2021 апрель
9 60 47 119. 2021 апрель
10 60 47.2 115. 2021 апрель
# … with 90,710 more rowsРассмотрим изменения основной величины wind_value для различных месяцев на гистограммах ниже.
Код
wind_2021 %>%
ggplot(aes(x = wind_value, fill = month)) +
geom_histogram(alpha = 0.9, binwidth = 15) +
labs(fill = "месяц",
x = "количество часов",
y = "количество наблюдений\n") +
scale_y_continuous(labels = function(x) format(x, big.mark = " ",
scientific = FALSE)) +
theme_bw() +
facet_wrap(~month, ncol = 2) +
theme(legend.position = "none",
plot.margin = margin(0.01, 0.01, 0.01, 0.01, "cm"),
legend.title = element_text(size = rel(1.1)),
legend.text = element_text(size = rel(1.1)))
Картографическая основа
Все карты (как векторные, так и растровые) и объекты на картах строятся с учетом картографической проекции, в нашем случае это проекция EPSG:4326. Кроме того, для геометрического упрощения границ и ускорения вычислений, заранее была выполнена генерализация картографической основы.
Интересующая нас область представляет собой прямоугольник, градусные меры широт которого варьируются от 42,75 до 60 и для долготы от 45 до 125,75.
Загрузим (упрощенные границы) Российской Федерации.
Код
Russia_map_path <- here(data_path, "boundary_Russia_regions_simplified.shp")
Russia_regions_boundary <- st_read(Russia_map_path)Reading layer `boundary_Russia_regions_simplified' from data source
`/Users/materov/@Мои доклады/@Quarto meetups/materials/report/data/boundary_Russia_regions_simplified.shp'
using driver `ESRI Shapefile'
Simple feature collection with 86 features and 1 field
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -180 ymin: 41.18587 xmax: 180 ymax: 82.05862
Geodetic CRS: WGS 84Код
Russia_regions_boundary <-
Russia_regions_boundary[order(Russia_regions_boundary$names_new), ]Изобразим полную карту РФ по регионам.
Код
ggplot() +
geom_sf(data = Russia_regions_boundary,
aes(fill = names_new),
alpha = 0.8) +
theme_void() +
theme(legend.position = "none") +
coord_sf(xlim = c(25, 180),
ylim = c(41.18587, 82.05862))
Ограничим карту регионов только на исследуемую прямоугольную область. Здесь цвет соответствует номеру региона.
Код
Russia_regions_crop <- st_crop(Russia_regions_boundary,
xmin = min(wind_2021$lon),
ymin = min(wind_2021$lat),
xmax = max(wind_2021$lon),
ymax = max(wind_2021$lat)
)
Russia_regions_crop <- Russia_regions_crop %>%
mutate(id = row_number())
ggplot() +
geom_sf(data = Russia_regions_crop,
aes(fill = id), alpha = 0.8) +
ggrepel::geom_label_repel(
data = Russia_regions_crop,
aes(label = id, geometry = geometry),
stat = "sf_coordinates",
min.segment.length = 0,
alpha = 0.7
) +
coord_sf(xlim = c(min(wind_2021$lon), max(wind_2021$lon)),
ylim = c(min(wind_2021$lat), max(wind_2021$lat)),
expand = FALSE) +
labs(x = "", y = "") +
theme(legend.position = "none") +
viridis::scale_fill_viridis(option = "turbo") +
theme(plot.margin = margin(0.01, 0.01, 0.01, 0.01, "cm"),
legend.title = element_text(size = rel(1.1)),
legend.text = element_text(size = rel(1.1)))
Замечание
Язык R позволяет строить график послойно, настраивать прозрачность слоев, добавлять различные элементы, и т.д.
Растровые карты
Для отображения данных можно воспользоваться растровыми картами с непрерывными палитрами. Ниже показан пример наблюдений июня 2021 года.
# функция, необходимая для переопределения
# средней точки цветовой палитры
1mid_rescaler <- function(mid = 0) {
function(x, to = c(0, 1), from = range(x, na.rm = TRUE)) {
scales::rescale_mid(x, to, from, mid)
}
}- 1
-
В коде выше функция
mid_rescaler– аргумент средней точки вscale_*(источник).
Код
library(paletteer)
ggplot() +
geom_raster(data = wind_2021 %>% dplyr::filter(month == "июнь"),
aes(x = lon, y = lat, fill = wind_value),
alpha = 0.9, interpolate = TRUE) +
geom_sf(data = Russia_regions_crop, color = "black",
alpha = 0.01,
linewidth = 0.4) +
ggrepel::geom_label_repel(
data = Russia_regions_crop,
aes(label = id, geometry = geometry),
stat = "sf_coordinates",
min.segment.length = 0,
alpha = 0.7
) +
scale_fill_paletteer_c("grDevices::RdYlBu",
na.value = NA,
rescaler = mid_rescaler(mid =
median(wind_2021$wind_value) + 100),
breaks = c(0, 100, 200, 300),
labels = c("0 ч", "100 ч", "200 ч", "300 ч"),
direction = -1) +
hrbrthemes::theme_ipsum() +
labs(x = "", y = "",
fill = "количество часов") +
coord_sf(xlim = c(min(wind_2021$lon), max(wind_2021$lon)),
ylim = c(min(wind_2021$lat), max(wind_2021$lat)),
expand = FALSE) +
guides(fill = guide_colorbar(title.position = 'top',
title.hjust = 0.5,
barwidth = unit(15, 'lines'),
barheight = unit(0.7, 'lines'))) +
theme(legend.position = "bottom",
plot.margin = margin(0.01, 0.01, 0.01, 0.01, "cm"),
legend.title = element_text(size = rel(1.1)),
legend.text = element_text(size = rel(1.1))
)
Карта высот над уровнем моря
Загрузим и визуализируем данные высот над уровнем моря, которые понадобятся нам для дальнейшего использования.
Код
# географическая проекция
projection <- "EPSG:4326"
# получение данных по высотам
elevation_df <- get_elev_raster(wind_2021 %>%
dplyr::select(lon, lat) %>%
rename(x = lon, y = lat) %>%
na.omit() %>%
as.data.frame(),
prj = projection,
# z = разрешение
z = 5)Из растровых данных легко сделать таблицу высот. Здесь \(x\) и \(y\) соответствуют координатам, а \(z\) – высотам над уровнем моря.
Код
elevation_tibble <- as.data.frame(elevation_df, xy = TRUE)
elevation_tibble %<>% as_tibble()
names(elevation_tibble)[3] <- "z"
elevation_tibble# A tibble: 4,165,928 × 3
x y z
<dbl> <dbl> <dbl>
1 45.0 61.6 160
2 45.0 61.6 167
3 45.1 61.6 173
4 45.1 61.6 172
5 45.1 61.6 161
6 45.1 61.6 147
7 45.1 61.6 138
8 45.2 61.6 143
9 45.2 61.6 155
10 45.2 61.6 161
# … with 4,165,918 more rowsПокажем растровое изображение высот над уровнем моря для исследуемой области.
Код
ggplot() +
geom_raster(data = elevation_tibble,
aes(x = x, y = y, fill = z)) +
geom_sf(data = Russia_regions_boundary, color = "black",
alpha = 0.1) +
ggrepel::geom_label_repel(
data = Russia_regions_crop,
aes(label = id, geometry = geometry),
stat = "sf_coordinates",
min.segment.length = 0,
alpha = 0.7
) +
scale_fill_etopo() +
hrbrthemes::theme_ipsum() +
labs(fill = "высота, м.:\n", x = "", y = "") +
theme(legend.key.height = unit(1, 'cm')) +
coord_sf(xlim = c(min(wind_2021$lon), max(wind_2021$lon)),
ylim = c(min(wind_2021$lat), max(wind_2021$lat)),
expand = FALSE) +
theme(plot.margin = margin(0.01, 0.01, 0.01, 0.01, "cm"),
legend.title = element_text(size = rel(1.1)),
legend.text = element_text(size = rel(1.1)))
Техническая информация
Здесь указана информация о сессии и загруженных библиотеках языка R.
Сноски
Желательно всегда представлять данные именно в таком виде.↩︎