ПЕРЕДМОВА

На самому початку, до вступного тексту, давайте в першу чергу подивимось параметри програмного забезпечення, яке буде використано нами упродовж нашого дослідження! Для цього використаємо базову функцію sessionInfo.

sessionInfo()
## R version 4.0.5 (2021-03-31)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Debian GNU/Linux 10 (buster)
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.8.0
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.8.0
## 
## locale:
##  [1] LC_CTYPE=uk_UA.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=uk_UA.UTF-8        LC_COLLATE=uk_UA.UTF-8    
##  [5] LC_MONETARY=uk_UA.UTF-8    LC_MESSAGES=uk_UA.UTF-8   
##  [7] LC_PAPER=uk_UA.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=uk_UA.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## loaded via a namespace (and not attached):
##  [1] digest_0.6.27     R6_2.5.0          jsonlite_1.7.2    magrittr_2.0.1   
##  [5] evaluate_0.14     rlang_0.4.11      stringi_1.5.3     jquerylib_0.1.4  
##  [9] bslib_0.2.4       rmarkdown_2.7     tools_4.0.5       stringr_1.4.0    
## [13] xfun_0.22         yaml_2.2.1        compiler_4.0.5    htmltools_0.5.1.1
## [17] knitr_1.33        sass_0.3.1

Як видно з виводу функції sessionInfo — ми працюємо на комп’ютері зі встановленою операційною системою Debian GNU/Linux 10 (buster) та використовуємо мову програмування R версії 4: R version 4.0.5 (2021-03-31).

Функція sessionInfo також виводить інформацію про локалізацію операційної системи: uk_UA.UTF-8. Окремими розділами функція sessionInfo надає інформацію про завантажені в поточній сесії пакети і про додаткові пакети, які не завантажено у поточній сесії, але які містяться у загальному просторі імен. Також функція інформує нас, що мова R “в курсі” про встановлені на комп’ютері бібліотеки BLAS та LAPACK — обчислювальні бібліотеки лінійної алгебри.

ВСТУП

Стрімкий розвиток мереж громадського екологічного моніторингу сприяє накопиченню величезних обсягів даних. Кожний пристрій щохвилинно надсилає до центральної бази даних результати вимірювань показників температури, атмосферного тиску, вологості, концентрацій забруднюючих речовин, радіаційного фону… Мета цієї статті — показати загальний підхід до отримання та впорядкування даних громадського екологічного моніторингу, що збираються на порталі Eco-City (https://eco-city.org.ua). У подальшому планується серія статей, в яких буде описано процес аналізу та візуалізації отриманих даних.

Кожен власник пристрою Eco-City має особистий кабінет на порталі https://eco-city.org.ua, де він може завантажити архівні дані свого пристрою за тривалий інтервал часу. Більше того, громадська ініціатива Eco-City збирає та зберігає дані інших мереж громадського екологічного моніторингу: SaveDnipro, AirPollution тощо. Відповідно — власники пристроїв цих виробників також можуть отримати архіви із даними спостережень своїх пристроїв із особистого кабінету на порталі Eco-City. У нашому випадку ми будемо використовувати дані пристрою ecorada1, який встановлено у м.Кривий Ріг на вул.Чеська. Цей пристрій — один із найперших в мережі громадського екологічного моніторингу Eco-City. Файли даних із результатами моніторингу, які використовуються упродовж цієї статті, можна отримати отут.

Використовувані програмні пакети R

Для обробки та аналізу даних ми будемо використовувати метапакет tidyverse та, частково, пакет reshape2. Для отримання даних метеорологічних спостережень буде використано пакет stationaRy.

Використання пакету stationaRy потребує стабільного інтернет-з’єднання!

Завантажуємо потрібні нам пакети:

library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──
## ✓ ggplot2 3.3.3     ✓ purrr   0.3.4
## ✓ tibble  3.1.1     ✓ dplyr   1.0.5
## ✓ tidyr   1.1.3     ✓ stringr 1.4.0
## ✓ readr   1.4.0     ✓ forcats 0.5.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library(stationaRy)

РОБОТА ІЗ ДАНИМИ

Особистий кабінет користувача на порталі Eco-City

Для можливості доступу до архівних результатів екологічного моніторингу на порталі Eco-City існують особисті кабінети користувачів.

Нас цікавлять дані пристрою ecorada1. Щоб їх отримати — ми повинні знайти цей пристрій у переліку, та зайти на його сторінку.

В особистому кабінеті на сторінці конкретного пристрою моніторингу ми можемо вибирати за який період треба завантажити дані: тиждень, місяць, квартал, рік. В усіх цих випадках буде братися інтервал ДО часу надання запиту: тиждень назад від надання запиту, місяць назад від надання запиту…, рік назад від надання запиту на завантаження архіву. Окрім вибору періоду завантаження архівних даних, в особистому кабінеті також міститься історія попередніх запитів.

Коли архів сформовано, користувачу слід оновити сторінку браузера, або вийти в попереднє меню (кнопка Назад до станцій) та знов повернутися до меню станції. Раджу в якості формату завантаження використовувати ZIP-файли, оскільки архівні таблиці первинних даних можуть досягати обсягу декількох сотень Мб.

В результаті виконання запиту ви отримаєте файл із подвійним розширенням CSV.ZIP, в якому міститься CSV-таблиця із результатами вимірювань та додатковою службовою інформацією. Назва файлу виглядає подібно до: 10 20201111140301 [2019-11-11|2020-11-11], де 10 - ID пристрою в центральній базі даних проекту EcoCity, 20201111140301 - дата та час генерації файлу архіву, [2019-11-11|2020-11-11] - період, за який містяться дані в отриманому архіві. В нашому випадку ми отримали та розпакували два архівні файли:

  • 10 20201111140301 [2019-11-11|2020-11-11].csv

  • 10 20210105180201 [2020-10-05|2021-01-05].csv

Перший файл містить дані пристрою моніторингу №10 (ecorada1) з 11-11-2019 до 11-11-2020 (річний інтервал). Другий файл містить дані за період від 05-10-2020 до 05-01-2021 (квартальний інтервал). Зверніть увагу, що дані в обох файлах перетинаються на часовій шкалі! Відповідно, в зведеній таблиці даних у нас будуть дубльовані значення, які треба вилучити.

Отримання даних метеорологічних спостережень

Окрім даних EcoCity, для нас надзвичайно корисними є дані про метеорологічні умови. В містах, які мають поблизу міжнародні аеропорти, можна отримати дані від цивільних авіаметеорологічних станцій. Ці дані із періодичністю 30 або 60 хвилин передаються в ефір та містять інформацію про погодні умови в районі аеропорту. Для інших населених пунктів також можливо отримати архіви метеорологічних спостережень на звичайних метеостанціях (щоправда - не для всіх наявних станцій!). Для отримання метеорологічних даних ми використаємо пакет stationaRy, який дозволяє розшукувати та вибірково отримувати дані метеорологічних спостережень з архіву NOAA - Національного управління океанічних і атмосферних досліджень США.

Єдина проблема, яка є на теперішній момент: пакет stationaRy не дозволяє отримати дані за 2021 рік. Проте, автор пакету працює над її виправленням. Ми ж будемо аналізувати результати вимірювань пристроєм ecorada1 упродовж 2020 року

Щоб отримати весь перелік метеорологічних станцій з архіву NOAA нам слід виконати команду get_station_metadata. Запишемо вивід цієї команди в об’єкт (таблицю) all.stations.

get_station_metadata() -> all.stations

Подивимось перші 10 рядків таблиці all.stations.

head(all.stations, 10)

Як бачимо, з таблиці all.stations, загалом в архіві NOAA міститься інформація з 29743 метеорологічних станцій.

glimpse(all.stations)
## Rows: 29,743
## Columns: 16
## $ id         <chr> "007018-99999", "007026-99999", "007070-99999", "008260-999…
## $ usaf       <chr> "007018", "007026", "007070", "008260", "008268", "008307",…
## $ wban       <chr> "99999", "99999", "99999", "99999", "99999", "99999", "9999…
## $ name       <chr> "WXPOD 7018", "WXPOD 7026", "WXPOD 7070", "WXPOD8270", "WXP…
## $ country    <chr> NA, "AF", "AF", NA, "AF", "AF", NA, NA, NA, NA, "NO", "NO",…
## $ state      <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ icao       <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, "ENRS", "ENJA", NA,…
## $ lat        <dbl> 0.000, 0.000, 0.000, 0.000, 32.950, 0.000, NA, NA, NA, NA, …
## $ lon        <dbl> 0.000, 0.000, 0.000, 0.000, 65.567, 0.000, NA, NA, NA, NA, …
## $ elev       <dbl> 7018.0, 7026.0, 7070.0, 0.0, 1156.7, 8318.0, NA, NA, NA, NA…
## $ begin_date <date> 2011-03-09, 2012-07-13, 2014-09-23, 2005-01-01, 2010-05-19…
## $ end_date   <date> 2013-07-30, 2017-08-22, 2015-09-26, 2010-09-20, 2012-03-23…
## $ begin_year <int> 2011, 2012, 2014, 2005, 2010, 2010, 2016, 2016, 2016, 2016,…
## $ end_year   <int> 2013, 2017, 2015, 2010, 2012, 2010, 2016, 2016, 2016, 2016,…
## $ tz_name    <chr> "Etc/GMT", "Etc/GMT", "Etc/GMT", "Etc/GMT", "Asia/Kabul", "…
## $ years      <list> <2011, 2013>, <2012, 2014, 2016, 2017>, <2014, 2015>, <200…

Щоб знайти метеорологічні станції України нам треба виконати серію команд. Запишемо результат в окрему таблицю ukr.stations.

get_station_metadata() %>%
  filter(country == "UP") -> ukr.stations

За результатами нашого пошуку видно, що в таблиці ukr.stations міститься інформація про 213 метеорологічних станцій (включно із тими, які знаходяться на тимчасово окупованих територіях).

Для нашого дослідження нам потрібно дані з метеорологічної станції, яка знаходиться в міжнародному аеропорту <<Кривий Ріг>>. Ця станція має в таблиці назву LOZUVATKA INTL та код ICAO UKDR.

Окрім основних метеорологічних даних, в архівах зберігаються багато додаткових даних щодо проведених метеорологічних спостережень. Різні метеорологічні станції мають різні набори додаткових даних. Детальний перелік та опис додаткових масивів даних викладено в документі від NOAA: https://www1.ncdc.noaa.gov/pub/data/ish/ish-format-document.pdf

Подивитися наявні набори додаткових даних можна з використанням допоміжної команди station_coverage. Збережемо вивід цієї команди в окрему таблицю KrR.coverage.

get_station_metadata() %>%
  dplyr::filter(name == "LOZUVATKA INTL") %>%
  dplyr::pull(id) %>%
  station_coverage(years = 2020) -> KrR.coverage

Загалом таблиця KrR.coverage містить інформацію про 87 додаткових наборів даних. Нас цікавить набір даних MA1, в якому містяться дані про атмосферний тиск. В процесі завантаження даних метеорологічних спостережень ми вкажемо що нам треба завантажити також і дані цього набору.

glimpse(KrR.coverage)
## Rows: 87
## Columns: 3
## $ id       <chr> "330003-99999", "330003-99999", "330003-99999", "330003-99999…
## $ category <chr> "AA1", "AB1", "AC1", "AD1", "AE1", "AG1", "AH1", "AI1", "AJ1"…
## $ count    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…

Тепер можна, власне, завантажити дані метеорологічних спостережень станції LOZUVATKA INTL за 2020 рік. При цьому ми примусово вимикаємо функцію make_hourly = FALSE щоб дані не усереднювались одногодинними інтервалами.

get_station_metadata() %>%
  dplyr::filter(name == "LOZUVATKA INTL") %>%
  dplyr::pull(id) %>%
  get_met_data(years = 2020,
               add_fields = c("MA1"),
               make_hourly = FALSE) -> KrR_metar.2020

Час метеорологічних спостережень в архівах NOAA записано в форматі UTC. Нам треба перевести час у київський часовий пояс.

KrR_metar.2020$time.Kyiv <- as.POSIXct(format(KrR_metar.2020$time,
      tz = "Europe/Kiev", usetz = TRUE))

Після отримання даних бажано зберегти їх у вигляді окремого архіву, щоб в подальшому не навантажувати ресурси NOAA повторними запитами.

saveRDS(KrR_metar.2020, "KrR_metar_2020.rds")

Обробка архівних даних Eco-City

Оскільки ми завантажили два файли з архівами екологічного моніторингу, ми повинні бути націленими на те, що доведеться автоматизована обробка цих файлів. Більше того — в реальному випадку ви можете мати й десяток таких файлів! Це означає, що ручне завантаження та ручний контроль буде вимагати багато часу та зусиль!

Щоб автоматизувати процес завантаження файлів можна створити в робочому просторі R окремий об’єкт, який буде містити назви всіх потрібних нам файлів. Для цього можна орієнтуватися на розширення файлів (у нашому випадку — CSV).

ecofiles <- dir(pattern = "*.csv")

Створений нами об’єкт ecofiles є звичайним текстовим вектором, який містить назви всіх CSV-файлів із робочої директорії. Ясна річ, що для того щоб такий підхід нормально працював — у робочій директорії не повинно бути інших CSV-файлів, окрім архівів екологічного моніторингу!

Тепер ми можемо завантажити всі наявні у створеному списку файли. Для цього ми використаємо команду read_csv пакету readr (входить до метапакету tidyverse). Після завантаження ми запустимо збирача сміття (команда gc — Garbage Collection). І, наприкінці процесу, бажано зберегти наше робоче середовище.

ecodata <- ecofiles %>%
  map(readr::read_csv, col_names = FALSE) %>%
  reduce(rbind)
## 
## ── Column specification ────────────────────────────────────────────────────────
## cols(
##   X1 = col_double(),
##   X2 = col_character(),
##   X3 = col_character(),
##   X4 = col_character(),
##   X5 = col_double(),
##   X6 = col_datetime(format = "")
## )
## 
## 
## ── Column specification ────────────────────────────────────────────────────────
## cols(
##   X1 = col_double(),
##   X2 = col_character(),
##   X3 = col_character(),
##   X4 = col_character(),
##   X5 = col_double(),
##   X6 = col_datetime(format = "")
## )
gc()
##            used  (Mb) gc trigger  (Mb) max used  (Mb)
## Ncells  1327171  70.9    2640128 141.0  2640128 141.0
## Vcells 27339001 208.6   69760890 532.3 65904729 502.9
save.image("EcoCity_RPubs.RData")

Отримана нами об’єднана таблиця містить 3155111 рядків та 6 стовпців. Давайте подивимось на її перші 10 рядків.

head(ecodata, 10)

Стовпець X1 — це ID пристрою в центральній базі даних порталу Eco-City. Стовпець X2 — внутрішня інформація про вимірювані параметри (для нас вона не є важливою). Стовпець X3 містить власне назву вимірюваного параметру. У стовпі X4 міститься інформація про вимірювані величини. А от у стовпці X5 міститься найцінніше — власне значення вимірювань! Стовпець X6 містить записи про дату-час. Зверніть увагу на те, що в процесі завантаження даних записи дати-часу було автоматично сконвертовано у формат POSIXct.

Отримана нами зведена таблиця представлена у так званому “довгому” форматі. Для більшості випадків аналізу даних нам треба її переформатувати у “широкий формат”. При цьому ми також вилучемо дубльовані записи, які обумовлені накладанням наших архівів на часовій шкалі. Дубльовані записи ми видаляємо простим обчисленням середнього арифметичного значення. Ясна річ, що можна було застосувати спеціальні функції для пошуку та видалення дубльованих записів. Для перетворення таблиці в довгу форму ми використаємо команду dcast пакету reshape2. Щоб не завантажувати цей пакет заради однократного виконання, ми просто вкажемо “повну” назву команди: reshape2::dcast. Після цих операцій ми також прибираємо сміття та зберігаємо файл робочого середовища.

ecodata %>%
  select(-c(X2, X4)) %>%
  reshape2::dcast(X1 + X6 ~ X3, value.var = "X5", fun.aggregate = mean, na.rm = TRUE) -> ecodata.wide
gc()
##            used  (Mb) gc trigger  (Mb) max used  (Mb)
## Ncells  1338218  71.5    6677116 356.6  8346395 445.8
## Vcells 31425940 239.8   69760890 532.3 69634189 531.3
save.image("EcoCity_RPubs.RData")

Тепер ми маємо таблицю ecodata.wide, яка представлена звичною для більшості людей “широкою” формою. Таблиця містить 580541 рядків та 7 стовпців. Давайте подивимось її перші 10 рядків.

head(ecodata.wide, 10)

Також можна подивитись загальні статистичні характеристики всіх стовпців таблиці.

summary(ecodata.wide)
##        X1           X6                         Humidity          PM10        
##  Min.   :10   Min.   :2019-11-11 00:00:00   Min.   :31.30   Min.   :   0.00  
##  1st Qu.:10   1st Qu.:2020-02-26 06:54:15   1st Qu.:96.90   1st Qu.:   7.17  
##  Median :10   Median :2020-06-09 08:45:50   Median :99.90   Median :  15.23  
##  Mean   :10   Mean   :2020-06-09 01:57:53   Mean   :93.17   Mean   :  32.60  
##  3rd Qu.:10   3rd Qu.:2020-09-20 23:53:19   3rd Qu.:99.90   3rd Qu.:  34.50  
##  Max.   :10   Max.   :2021-01-05 18:01:32   Max.   :99.90   Max.   :1999.90  
##                                             NA's   :173     NA's   :42       
##      PM2.5         Temperature      WiFiSignal    
##  Min.   :  0.00   Min.   :-9.90   Min.   :-87.00  
##  1st Qu.:  3.33   1st Qu.: 3.60   1st Qu.:-70.00  
##  Median :  7.70   Median :11.20   Median :-65.00  
##  Mean   : 14.14   Mean   :12.63   Mean   :-66.14  
##  3rd Qu.: 18.42   3rd Qu.:20.70   3rd Qu.:-62.00  
##  Max.   :999.90   Max.   :40.50   Max.   : 31.00  
##  NA's   :43       NA's   :173     NA's   :1

Ми бачимо, що часовий інтервал вимірювань охоплює проміжок часу від 2019-11-11 до 2021-01-05 18:01:32. При цьому ми маємо 173 порожні значення (NA) відносної вологості повітря (Humidity), 42 порожні значення концентрацій пилу PM10, 43 порожні значення концентрації пилу PM2.5, 173 порожні значення показників температури повітря (Temperature) та 1 порожнє значення потужності WiFi-сигналу. При загальній довжині таблиці 580541 рядків — зовсім мізерні показники втрачених даних!