Фон

В Зимбабве субтропический высокогорный климат с сезоном дождей с ноября по март и сухим сезоном с мая по сентябрь. На протяжении многих лет Зимбабве сталкивалась с различными проблемами, связанными с климатом, такими как засухи, аномальная жара и наводнения, которые оказали значительное влияние на экономику и население страны.

Недавно в марте 2019 года на Зимбабве обрушился циклон Идай, который стоил стране 622 миллиона долларов США и затронул 270 000 человек. Ущерб был нанесен инфраструктуре, домам, посевам, домашнему скоту и другим ресурсам. Всемирный банк сообщил, что для лучшего восстановления предполагаемые затраты составят 1,1 миллиарда долларов США.

Теперь вопрос, который волнует многих политиков, страховых компаний и ученых, заключается в том, как мы можем проанализировать климатические данные Зимбабве, чтобы лучше понять эти проблемы и информировать о соответствующей политике и мерах по адаптации?

Что ж, мы можем ответить на этот вопрос, проанализировав различные климатические показатели, такие как температура, осадки и экстремальные погодные явления, за определенный период, скажем, 30 лет. Это позволит нам определить тенденции и изменения климата и их потенциальное воздействие на различные секторы экономики, такие как сельское хозяйство, энергетика и водные ресурсы. Мы также можем использовать климатические модели и сценарии для прогнозирования будущих климатических условий и предоставления информации для долгосрочного планирования и принятия решений. Понимая динамику климата в Зимбабве, мы можем разработать основанные на фактических данных стратегии для смягчения климатических рисков и повышения устойчивости, а также для содействия устойчивому развитию.

Цель

Цель этого блога — научить читателя тому, как моделировать климатические показатели Зимбабве на уровне провинций с использованием некоторых моделей машинного обучения.

Методология

Для достижения нашей цели мы будем следовать конвейеру рабочего процесса.

Пакет target в R имеет много преимуществ, поскольку создает рабочий процесс. То, что вы видите выше, — это взаимосвязанные функции, которые выполняют задачи. Эти задачи будут рассмотрены. С помощью этого рабочего процесса 1. вы можете увидеть, какой код не работает, и можете быстро решить проблему, 2. вы сократите время, затрачиваемое на выполнение кода, поскольку этот рабочий процесс будет повторно запускать только измененный код.

Сбор данных

Итак, сначала мы получим данные из двух API. Что такое API? API расшифровывается как «интерфейс прикладного программирования». Это набор протоколов, подпрограмм и инструментов, используемых разработчиками программного обеспечения для создания программных приложений.

Сначала мы получим данные о долготе и широте с этого сайта: https://nominatim.openstreetmap.org/. Эта функция принимает провинцию/место (Хараре, Булавайо и т. д.) как x и приклеивает/прикрепляет это место к URL-адресу API, затем взаимодействует с файлом JSON и затем выводит таблицу долготы и широты. Обратите внимание: для пробелов, например (Mashonaland East), я удалил пробел между ними и заменил их на +. Однако API может изменить свои требования к вводу, так что помните об этом.

get_location <- function(x){
  
  url <- glue::glue("https://nominatim.openstreetmap.org/search?q={x}&format=json")
  
  # Retrieve search results as JSON
  results <- jsonlite::fromJSON(url)
  
  # Extract latitude and longitude from first result
  lat <- results$lat[1]
  lon <- results$lon[1]
  
  # save as a tibble
  tibble(
    lat, lon
  ) 
  
}

Что касается широты и долготы, наша следующая задача — получить данные с 1992 по 2022 год. Для этого упражнения мы возьмем наши данные с этого сайта https://archive-api.open-meteo.com/v1/archive?. Снова склеиваем наши значения широты и долготы и получаем следующие показатели: 1. Максимальная температура, 2. Средняя температура, 3. Минимальная температура, 4. Сумма осадков, 5. Количество осадков в часах, 6. Сумма дождей. Мы снова сохраняем это как таблицу.

get_data <- function(lat, lon){
  
  start = "1992-01-01"
  end = "2022-12-31"
  
  # Website where we getting the data from
  url = glue::glue("https://archive-api.open-meteo.com/v1/archive?latitude={lat}&longitude={lon}&start_date={start}&end_date={end}&daily=temperature_2m_max,temperature_2m_min,temperature_2m_mean,precipitation_sum,rain_sum,precipitation_hours&timezone=Africa%2FCairo")
  
  # Get the data from the JSON file
  response <- jsonlite::fromJSON(url) 
  
  # get the variables
  date <- response$daily$time
  max_temperature <- response$daily$temperature_2m_max
  mean_temperature <- response$daily$temperature_2m_mean
  min_temperature <- response$daily$temperature_2m_min
  precipitation_sum <- response$daily$precipitation_sum
  rain_sum <- response$daily$rain_sum
  precipitation_hours <- response$daily$precipitation_hours
  
  # Save as a tibble
  tibble(
    date, max_temperature, mean_temperature, min_temperature, precipitation_sum,
    rain_sum, precipitation_hours
  )
  
}

Чтобы получить наши данные, нам нужно использовать функцию pmap в пакете purrr, который также находится в пакете tidyverse. Обратите внимание, что ниже используются функции get_location и get_data. Мы сохраняем наши данные как файл rds, почему? ответ почему нет?

get_all_df <- function(province){
  df <- province %>% 
    as_tibble() %>% 
    mutate(loc_stat = map(.x = provinces, 
                          .f = ~get_location(.x))) %>% 
    unnest(loc_stat) %>% 
    mutate(data = pmap(list(lat, lon), get_data)) %>% 
    unnest(data) %>% 
    mutate(date = as.Date(date)) %>% 
    mutate(day = wday(date, label = T),
           year = year(date),
           week = week(date))
  
  # Save data
  saveRDS(df, "output/climate_data.rds")
  
  tibble(df)
}

Иногда здорово хранить данные в виде базы данных. Отсюда любой, у кого есть навыки работы с SQL, может играть с данными. Делаем это следующим образом:

store_data_db <- function(data){
  # create a connection to the database
  con <- dbConnect(RSQLite::SQLite(), 
                   "output/climate_data.db")
  
  # create a table in the database and insert the data
  dbWriteTable(con, 
               "zimbabwe_table", 
               data, 
               overwrite = TRUE)
  
  
  # close the connection
  dbDisconnect(con)
}

Анализ данных

Так как же выглядят данные? мы можем выбрать, скажем, Булавайо, город королев и королей, и индикатор средней температуры, чтобы посмотреть на данные. Вы можете просто удалить prov/var и заменить фактический код. Этот код был подготовлен для среды Shiny Dashboard.

subset_data <- function(data){
  # Province
  prov <- "Bulawayo"
  
  # Variable
  var <- "mean_temperature"
  
  # Data
  subset_climate <- data %>% 
    filter(provinces  == prov) %>% 
    select(date, var)
  
}

trend_plot <- function(data){
  
  # Plot the data
  data %>% 
    mutate(date = as.Date(date)) %>% 
    plot_time_series(date, mean_temperature) 
}

С этим кодом мы получаем красивый сюжет. Мы можем наблюдать тенденцию и то, что в среднем есть жаркие и холодные дни.

Чтобы оценить точность наших моделей, нам нужно разделить данные на обучающую и тестовую выборки. С помощью тестового набора мы можем увидеть, как наши обученные модели работают с данными, которых они не видели. Мы будем использовать последние 4 месяца в качестве тестового набора и начертим наш план.

split_data <- function(data){
  split <- time_series_split(
    data,
    assess = "4 months",
    cumulative = TRUE
  )
}


trend_split_plot <- function(data){
  
  data %>% 
    tk_time_series_cv_plan() %>% 
    plot_time_series_cv_plan(date, mean_temperature)
  
  
  
}

План — это небольшой кусок по сравнению с остальными данными

Теперь, что касается моделей, мы будем обучать 1. Модель авторегрессионного интегрированного скользящего среднего (ARIMA), 2. Модель Пророка и, наконец, 3. Модель GLMNET. Последняя модель получит месяц, неделю и день для моделирования данных.

modeling <- function(data){
  
  ##Auto Arima
  
  model_arima <- arima_reg() %>% 
    set_engine("auto_arima") %>% 
    fit(mean_temperature ~ date, training(data))
  
  
  ## Prophet
  
  model_prophet <- prophet_reg() %>% 
    set_engine("prophet") %>% 
    fit(mean_temperature ~ date, training(data))
  
  ## GLMNET
  model_glmnet<- linear_reg(penalty = 0.01) %>% 
    set_engine("glmnet") %>% 
    fit(mean_temperature ~ wday(date, label = T)
        + month(date, label = T)
        + as.numeric(date), 
        training(data))
  
  # Table
  model_tbl <- modeltime_table(
    model_arima,
    model_prophet,
    model_glmnet
  )
  
  
}

Получаем следующие результаты:

Судя по среднеквадратической ошибке (RMSE), мы можем заметить, что модель пророка превосходит модель Seasonal ARIMA (SARIMA) и модели GLMNET.

Мы можем получить тестовые изображения, используя следующую функцию

get_test_viz <- function(data, split_data, one_province){
  
  data %>% 
    modeltime_forecast(
      new_data = testing(split_data),
      actual_data = one_province
    ) %>% 
    plot_modeltime_forecast()
}

Как и ожидалось, мы получаем следующие изображения, показывающие, как модель пророка пытается уловить тенденцию.

Далее мы делаем некоторые прогнозы. С помощью этой функции мы прогнозируем будущее на 4 месяца и получаем график

get_forecast <- function(data, one_province){
  
  future_forecast_tbl <- data %>% 
    modeltime_refit(one_province) %>% 
    modeltime_forecast(
      h = "4 months",
      actual_data = one_province
    )

Наконец, в качестве бонуса, как мы можем сделать такую ​​картинку?

Мы очищаем знак +, который мы включили ранее, чтобы он читался как Mat. Юг, а не Мат + Юг, это упростит соединение с шейп-файлом. Итак, где я взял шейп-файл? Чтобы узнать больше об этом, я создам еще один блог о создании карт и о том, где взять шейп-файлы.

Следующий код выбирает 4 января 2012 года и отображает среднюю температуру. Вы можете выбрать любую дату/подмножество диапазона/любой индикатор и посмотреть, как все пойдет.

get_map_data <- function(data){
  # Select variable
  vars = "mean_temperature"
  
  # Select the date
  date_select = "2012-01-04"
  
  climate_data <- data %>%
    mutate(provinces = str_replace(provinces, "\\+", " ")) %>% 
    filter(date == date_select)  
  
  # Zim Shape file
  zim_sh <- st_read("data/zimbabwe_shape_files/ZWE_adm1.shp")
  
  # Combine the data and shape file
  zim_sh <- zim_sh %>%
    left_join(., climate_data, by = c("NAME_1" = "provinces"))
  
}


get_map <- function(data){
  
  
  # Select the palette
  palette = "YlGn"
  
  # Create the plot
  zim_map <- data %>%
    ggplot() +
    geom_sf(aes(fill = mean_temperature)) +
    scale_fill_distiller(palette = palette, trans = "reverse") +
    theme_void()
  
  # Plotly map
  plotly::ggplotly(zim_map)
  
}

Заключение

Нам удалось сделать следующее:

  1. Получить данные из источников API
  2. Использовали конвейер рабочего процесса, чтобы сделать нашу историю правильной.
  3. Смоделируйте данные Булавайо, используя ARIMA, Prophet и GLMNET, из которых модель Prophet превзошла другие модели.
  4. Мы сделали график ландшафта Зимбабве, показывающий, что Мэт-Норт и Масвинго имеют самые высокие средние температуры.

Вооружившись такими моделями, мы можем заранее знать о резких климатических условиях впереди и мы можем быстро предпринять необходимые шаги для управления надвигающимися опасностями.

Код можно найти здесь: Github