Работа с временными рядами или пространственными данными подразумевает наличие “генеральной линии” вдоль которой происходит изменение в данных. В случае временного ряда это будет координата времени, в случае пространственных данных – какое-либо пространственное направление. Но в любом наборе данных, независимо от масштаба – будь то микроскопический или космический, мы имеем дискретность и пропуски в данных. Дискретность объясняется физической невозможностью исследовать непрерывно и с высокой детальностью интересующие нас величины. Пропуски в данных обусловлены всякими техническими неполадками – сбоями в работе оборудования, утратой документов, ошибками персонала при сохранении или копировании и т.д.
Иногда приходится объединять разные наборы данных. Проблем не возникает, если они сведены в едином масштабе (разрешении) и с одинаковой детальностью. Но что делать – если данные из абсолютно разных областей исследований и сами исследования проводились в разных масштабах? В этой статье мы рассмотрим именно такой вариант – сведение воедино разных наборов данных, записанных в разных масштабах. Для решения этой задачи мы воспользуемся функциями пакетов zoo и dplyr. Для анализа также применим возможности dplyr, а для отображения результатов и получения красивой графики – используем пакет ggplot2. При решении всех задач мы будем использовать туннелирование, которое позволяет не плодить большое количество промежуточных таблиц.
Начнём!
У нас есть данные мониторинга качества воздуха и данные метеорологических наблюдений в городе Кривой Рог. Нам нужно проанализировать зависимости между метеоусловиями и загрязнением воздуха. В наборе данных мониторинга качества воздуха присутствуют названия постов мониторинга (nazva), дата-время измерений (time.posCT.Kiev) и значения концентраций пыли (dust). В наборе метеорологических данных присутствуют дата-время измерений (time.posCT.Kiev), значения скорости ветра (ws) и значения направления ветра (wd). Направление ветра указано в градусах и “откуда дует”.
Мониторинг качества воздуха производится автоматизированными постами с интервалом 20 минут: в 00, 20 и 40 минут каждого часа. Метеорологические наблюдения производятся с интервалом 30 минут: в 00 и 30 минут каждого часа. Что делать? – выбрасывать часть данных, смещать интервалы, или растягивать интервалы??? Поскольку основным объектом исследования является загрязнение воздуха, то именно к этим данным мы будем “присоединять” все остальные наборы данных. И именно данные о загрязнении нам необходимо сохранить, по возможности, в неизменном виде! Поэтому все преобразования мы будем проводить с данными метеорологических наблюдений.
“Нулевым этапом” нашей работы будет подготовка данных: загрузка и приведение времени к формату POSIXct. Именно с этим форматом без проблем работает пакет dplyr.
atmo2 <- read.csv("atmo2.csv", header = T, sep = ";", dec = ".")# загружаем таблицу с данными мониторинга качества воздуха
KRmeteo <- read.csv("KRmeteo.csv", header = T, sep = ";", dec = ".")# загружаем таблицу с данными метеорологических наблюдений
atmo2$time.posCT.Kiev <- as.POSIXct(as.character(atmo2$time.posCT.Kiev), format = "%Y-%m-%d %H:%M:%S", tz = "Europe/Kiev")# приводим время мониторинга в формат POSIXct
KRmeteo$time.posCT.Kiev <- as.POSIXct(as.character(KRmeteo$time.posCT.Kiev), format = "%Y-%m-%d %H:%M:%S", tz = "Europe/Kiev")# приводим время метеонаблюдений в формат POSIXct
Теперь займёмся собственно работой!
На первом этапе нам нужно сделать так, чтобы временные интервалы метеорологических наблюдений совпадали с интервалами мониторинга качества воздуха. Самый простой вариант – провести интерполяцию метеоданных и высчитать промежуточные значения в интервалах 10, 20, 40 и 50 минут каждого часа (не забываем, что первичные данные содержат измерения в 00 и 30 минут!). Для интерполяции мы воспользуемся функцией na.approx из пакета zoo и проведём линейную интерполяцию:
library(zoo)#загружаем пакет zoo KRmeteo$wd[is.na(KRmeteo$wd)] <- -9999# устанавливаем для отсутствующих значений "минус почти десять тысяч", зачем - потом объясню seq.time <- zoo(order.by=(as.POSIXct(seq(min(KRmeteo$time.posCT.Kiev), max(KRmeteo$time.posCT.Kiev), by=600))))# создаём упорядоченный временной ряд с шагом 600 секунд mer.meteo <- merge(zoo(x = KRmeteo[c(11:15, 17)], order.by = unique(KRmeteo$time.posCT.Kiev)), seq.time)# объединяем первичную таблицу с новым временным рядом approx.meteo <- round(na.approx(mer.meteo),0)# проводим интерполяцию временного ряда и округляем полученные значения до целого
Для преобразования временного ряда из служебного формата пакета zoo в простой data.frame – нам понадобится специальная функция:
zooToDf <- function(z) {# создаём функцию для преобразования временного ряда в таблицу данных
df <- as.data.frame(z)
df$Date <- time(z) # создаём поле даты-времени
rownames(df) <- NULL# без заголовков строк
df <- df[,c(ncol(df), 1:(ncol(df)-1))]# упорядочиваем поля и ставим поле даты-веремени на первое место
return(df)
}
approx.meteo.df <- zooToDf(approx.meteo)# преобразуем временной ряд в таблицу данных
Теперь у нас есть таблица данных метеонаблюдений с интервалом в 10 минут. Именно эту таблицу мы будем присоединять к таблице мониторинга качества воздуха. Поскольку в реальности оба набора данных довольно большие – функция merge может не справиться с такими объёмами. Для объединения таблиц мы воспользуемся функцией left_join пакета dplyr. Но вначале мы пересчитаем направление ветра с градусов на румбы. Воспользуемся 16-румбовой схемой. И тут объясним – зачем мы приписывали отрицательное значение пустым записям. Дело в том, что помимо ветра определённого направления возможны ещё два варианта метеоусловий: слабый ветер переменного направления и полный штиль. Для слабого ветра не будет значения в поле “направление”, но будет какое-то число в поле “скорость”. Для полного штиля в поле “скорость” будет “0”, а в поле “направление” – будет пустое значение. При интерполяции мы отсекаем все возможные случаи переменного ветра и штиля между известными значениями. Значение направления ветра в румбах мы создадим в виде упорядоченного фактора и вручную зададим порядок значений (с севера по часовой стрелке). В самый хвост фактора мы поместим значения “переменный ветер” (VRB) и “полный штиль” (CALM). Упорядочение снимет в дальнейшем головную боль при построении графиков. При пересчёте воспользуемся туннелированием и функцией mutate пакета dplyr.
library(dplyr)# загружаем пакет dplyr
tbl_df(approx.meteo.df) %>% mutate(DD.meteo.from=ordered(ifelse((wd >= 168.75 & wd < 191.25 & ws > 0),"S", ifelse((wd >= 191.25 & wd < 213.75 & ws > 0),"SSW", ifelse((wd >= 213.75 & wd < 236.25 & ws > 0),"SW", ifelse((wd >= 236.25 & wd < 258.75 & ws > 0),"WSW", ifelse((wd >= 258.75 & wd < 281.25 & ws > 0),"W", ifelse((wd >= 281.25 & wd < 303.75 & ws > 0),"WNW", ifelse((wd >= 303.75 & wd < 326.25 & ws > 0),"NW", ifelse((wd >= 326.25 & wd < 348.75 & ws > 0),"NNW", ifelse((wd >= 348.75 & wd <= 360),"N", ifelse((wd >= 0 & wd < 11.25 & ws > 0),"N", ifelse((wd >= 11.25 & wd < 33.75 & ws > 0),"NNE", ifelse((wd >= 33.75 & wd < 56.25 & ws > 0),"NE", ifelse((wd >= 56.25 & wd < 78.75 & ws > 0),"ENE", ifelse((wd >= 78.75 & wd < 101.25 & ws > 0),"E", ifelse((wd >= 101.25 & wd < 123.75 & ws > 0),"ESE", ifelse((wd >= 123.75 & wd < 146.25 & ws > 0),"SE", ifelse((wd >= 146.25 & wd < 168.75 & ws > 0),"SSE", ifelse((wd < 0 & ws > 0), "VRB", "CALM")))))))))))))))))), levels=c("N", "NNE", "NE", "ENE", "E", "ESE", "SE", "SSE", "S", "SSW", "SW", "WSW", "W", "WNW", "NW", "NNW", "VRB", "CALM"))) -> approx.meteo.df# производим все вычисления и записываем результат обратно в таблицу approx.meteo.df
Теперь объединим обе таблицы:
atmo2 <- left_join(atmo2, approx.meteo.df)# объединяем две таблицы данных
Функция left_join сама выберет одинаковые поля, по которым и произведёт объединение таблиц. В нашем случае это будут поля time.posCT.Kiev в обеих таблицах. Других совпадающих названий в таблицах быть не должно! Иначе на выходе вы получите сложную таблицу с агрегированными по нескольким ключам данными.
Теперь давайте сделаем с нашими данными что-нибудь умное! Например – построим диаграмму Тьюки, на которой будет отображаться изменчивость концентраций пыли при разных направлениях ветра.
png("WS_dust.png", width = 420, height = 297, units = "mm", res = 150)# сразу запускаем экспорт во внешний файл
tbl_df(atmo2) %>%# запускаем в туннель таблицу с данными
select(nazva, DD.meteo.from, dust) %>%# отбираем нужные для построения диаграммы поля
ggplot(aes(y = dust, x = DD.meteo.from, fill = nazva))+# задаём общие параметры диаграммы
geom_boxplot()+# задаём тип диаграммы
stat_summary(fun.y=median, geom="line", aes(group=2), colour = "red", size = 2)+# строим линию, объединяющую медианы на диаграмме
scale_y_continuous(limits=c(0.0, 0.75))+# устанавливаем границу значений по оси Y (не отображаем сильно высокие значения концентраций)
guides(fill = FALSE)+# убираем отображение легенды к диаграмме (она не нужна - названия постов отображаются на самой диаграмме)
facet_grid(nazva~.)# разделяем диаграмму по названиям постов мониторинга качества воздуха
dev.off()# закрываем графический интерфейс
Зависимость концентраций пыли в воздухе от направления ветра
И вуаля!!! Теперь у нас есть сводный график по всем трём постам, на котором отображена изменчивость концентраций пыли в воздухе в зависимости от направления ветра (“откуда дует”). Конечно, этот график можно ещё допилить – добавить подписи, шрифты и прочие рюшечки… Но это уже совсем другая история… 🙂
PS Примеры оформления графиков Вы можете посмотреть в предыдущих статьях на этом сайте. Например, в этой статье.