# 1: Load the packages 

library(readr)
library(dplyr)
library(ggplot2)
library(tidyr)
library(knitr)
# 2: Import the meteorology data file
# I chose the Seyakha station for this project.

weather_raw <- read_fwf(
  file = "uni.seyakha.20967.dat",
  col_positions = fwf_empty("uni.seyakha.20967.dat")
)
## Rows: 480 Columns: 22
## ── Column specification ────────────────────────────────────────────────────────
## 
## chr  (2): X9, X22
## dbl (20): X1, X2, X3, X4, X5, X6, X7, X8, X10, X11, X12, X13, X14, X15, X16,...
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
# Look at the first few rows to check the data
head(weather_raw)
## # A tibble: 6 × 22
##      X1    X2    X3    X4    X5    X6    X7    X8 X9      X10   X11   X12   X13
##   <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 20967  1961     1    -1    -1     1  70.2  72.6 -25.4 1015. 1000. 1000.   6  
## 2 20967  1961     2    -1    -1     1  70.2  72.6 -24.1 1007. 1000. 1000.   6.8
## 3 20967  1961     3    -1    -1     1  70.2  72.6 -12.6 1004. 1000. 1000.   7.6
## 4 20967  1961     4    -1    -1     1  70.2  72.6 -13.3 1006  1000. 1000.   6.6
## 5 20967  1961     5    -1    -1     1  70.2  72.6 -8.0  1010. 1000. 1000.   9  
## 6 20967  1961     6    -1    -1     1  70.2  72.6 .6    1011. 1000. 1000.   8.7
## # ℹ 9 more variables: X14 <dbl>, X15 <dbl>, X16 <dbl>, X17 <dbl>, X18 <dbl>,
## #   X19 <dbl>, X20 <dbl>, X21 <dbl>, X22 <chr>
# 3: Check the column names and structure
# This helps me figure out which columns contain year, month, latitude, longitude, and air temperature.

names(weather_raw)
##  [1] "X1"  "X2"  "X3"  "X4"  "X5"  "X6"  "X7"  "X8"  "X9"  "X10" "X11" "X12"
## [13] "X13" "X14" "X15" "X16" "X17" "X18" "X19" "X20" "X21" "X22"
head(weather_raw)
## # A tibble: 6 × 22
##      X1    X2    X3    X4    X5    X6    X7    X8 X9      X10   X11   X12   X13
##   <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 20967  1961     1    -1    -1     1  70.2  72.6 -25.4 1015. 1000. 1000.   6  
## 2 20967  1961     2    -1    -1     1  70.2  72.6 -24.1 1007. 1000. 1000.   6.8
## 3 20967  1961     3    -1    -1     1  70.2  72.6 -12.6 1004. 1000. 1000.   7.6
## 4 20967  1961     4    -1    -1     1  70.2  72.6 -13.3 1006  1000. 1000.   6.6
## 5 20967  1961     5    -1    -1     1  70.2  72.6 -8.0  1010. 1000. 1000.   9  
## 6 20967  1961     6    -1    -1     1  70.2  72.6 .6    1011. 1000. 1000.   8.7
## # ℹ 9 more variables: X14 <dbl>, X15 <dbl>, X16 <dbl>, X17 <dbl>, X18 <dbl>,
## #   X19 <dbl>, X20 <dbl>, X21 <dbl>, X22 <chr>
# 4: Rename the important columns
# These are the columns I need for the project:
# station, year, month, day, latitude, longitude, and temperature.

weather <- weather_raw %>%
  rename(
    station = X1,
    year = X2,
    month = X3,
    day = X4,
    latitude = X7,
    longitude = X8,
    temperature_c = X9
  ) %>%
  mutate(
    temperature_c = as.numeric(temperature_c)
  )

head(weather)
## # A tibble: 6 × 22
##   station  year month   day    X5    X6 latitude longitude temperature_c   X10
##     <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>    <dbl>     <dbl>         <dbl> <dbl>
## 1   20967  1961     1    -1    -1     1     70.2      72.6         -25.4 1015.
## 2   20967  1961     2    -1    -1     1     70.2      72.6         -24.1 1007.
## 3   20967  1961     3    -1    -1     1     70.2      72.6         -12.6 1004.
## 4   20967  1961     4    -1    -1     1     70.2      72.6         -13.3 1006 
## 5   20967  1961     5    -1    -1     1     70.2      72.6          -8   1010.
## 6   20967  1961     6    -1    -1     1     70.2      72.6           0.6 1011.
## # ℹ 12 more variables: X11 <dbl>, X12 <dbl>, X13 <dbl>, X14 <dbl>, X15 <dbl>,
## #   X16 <dbl>, X17 <dbl>, X18 <dbl>, X19 <dbl>, X20 <dbl>, X21 <dbl>, X22 <chr>
# 5: Clean the temperature data
# Missing temperature values are coded as 99.99 or 999.99 so I removed them from the data.

weather_clean <- weather %>%
  filter(
    temperature_c != 99.99,
    temperature_c != 999.99,
    !is.na(temperature_c)
  )

head(weather_clean)
## # A tibble: 6 × 22
##   station  year month   day    X5    X6 latitude longitude temperature_c   X10
##     <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>    <dbl>     <dbl>         <dbl> <dbl>
## 1   20967  1961     1    -1    -1     1     70.2      72.6         -25.4 1015.
## 2   20967  1961     2    -1    -1     1     70.2      72.6         -24.1 1007.
## 3   20967  1961     3    -1    -1     1     70.2      72.6         -12.6 1004.
## 4   20967  1961     4    -1    -1     1     70.2      72.6         -13.3 1006 
## 5   20967  1961     5    -1    -1     1     70.2      72.6          -8   1010.
## 6   20967  1961     6    -1    -1     1     70.2      72.6           0.6 1011.
## # ℹ 12 more variables: X11 <dbl>, X12 <dbl>, X13 <dbl>, X14 <dbl>, X15 <dbl>,
## #   X16 <dbl>, X17 <dbl>, X18 <dbl>, X19 <dbl>, X20 <dbl>, X21 <dbl>, X22 <chr>
# 6: Create the two time periods
# I am comparing 1961-1980 with 1981-2000.

weather_clean <- weather_clean %>%
  mutate(
    period = case_when(
      year >= 1961 & year <= 1980 ~ "1961-1980",
      year >= 1981 & year <= 2000 ~ "1981-2000",
      TRUE ~ "Other"
    )
  ) %>%
  filter(period != "Other")

table(weather_clean$period)
## 
## 1961-1980 1981-2000 
##       240       240
# 7: Calculate monthly summary statistics
# This gives the mean, median, and standard deviation of temperature for each month in each time period.

monthly_summary <- weather_clean %>%
  group_by(period, month) %>%
  summarise(
    mean_temp = mean(temperature_c, na.rm = TRUE),
    median_temp = median(temperature_c, na.rm = TRUE),
    sd_temp = sd(temperature_c, na.rm = TRUE),
    n = n(),
    .groups = "drop"
  )

monthly_summary
## # A tibble: 24 × 6
##    period    month mean_temp median_temp sd_temp     n
##    <chr>     <dbl>     <dbl>       <dbl>   <dbl> <int>
##  1 1961-1980     1    -26.7       -26.0     2.76    20
##  2 1961-1980     2    -26.6       -26.2     5.44    20
##  3 1961-1980     3    -22.4       -22.0     4.81    20
##  4 1961-1980     4    -16.1       -16.6     3.67    20
##  5 1961-1980     5     -7.43       -7.3     1.78    20
##  6 1961-1980     6      1.04        1       1.23    20
##  7 1961-1980     7      7.27        7.4     1.65    20
##  8 1961-1980     8      7.44        8       1.48    20
##  9 1961-1980     9      3.29        3.05    1.29    20
## 10 1961-1980    10     -6.69       -6.85    2.95    20
## # ℹ 14 more rows
# 8: Make a table for the monthly summary statistics
# This table shows the monthly mean and median temperature for 1961-2000, 1961-1980, and 1981-2000.

full_period_summary <- weather_clean %>%
  group_by(month) %>%
  summarise(
    mean_1961_2000 = mean(temperature_c, na.rm = TRUE),
    median_1961_2000 = median(temperature_c, na.rm = TRUE),
    .groups = "drop"
  )

period_summary_wide <- monthly_summary %>%
  select(period, month, mean_temp, median_temp) %>%
  mutate(
    period = recode(
      period,
      "1961-1980" = "early",
      "1981-2000" = "later"
    )
  ) %>%
  pivot_wider(
    names_from = period,
    values_from = c(mean_temp, median_temp)
  )

monthly_summary_table <- full_period_summary %>%
  left_join(period_summary_wide, by = "month") %>%
  mutate(
    month_name = factor(
      month,
      levels = 1:12,
      labels = c("Jan", "Feb", "Mar", "Apr", "May", "Jun",
                 "Jul", "Aug", "Sep", "Oct", "Nov", "Dec")
    ),
    mean_1961_2000 = round(mean_1961_2000, 2),
    median_1961_2000 = round(median_1961_2000, 2),
    mean_temp_early = round(mean_temp_early, 2),
    median_temp_early = round(median_temp_early, 2),
    mean_temp_later = round(mean_temp_later, 2),
    median_temp_later = round(median_temp_later, 2)
  ) %>%
  select(
    month_name,
    mean_1961_2000,
    median_1961_2000,
    mean_temp_early,
    median_temp_early,
    mean_temp_later,
    median_temp_later
  )

kable(
  monthly_summary_table,
  col.names = c(
    "Month",
    "Mean 1961-2000",
    "Median 1961-2000",
    "Mean 1961-1980",
    "Median 1961-1980",
    "Mean 1981-2000",
    "Median 1981-2000"
  ),
  caption = "Table 1. Monthly mean and median air temperature statistics for 1961-2000, 1961-1980, and 1981-2000."
)
Table 1. Monthly mean and median air temperature statistics for 1961-2000, 1961-1980, and 1981-2000.
Month Mean 1961-2000 Median 1961-2000 Mean 1961-1980 Median 1961-1980 Mean 1981-2000 Median 1981-2000
Jan -25.80 -25.85 -26.71 -25.95 -24.90 -25.65
Feb -25.92 -25.30 -26.57 -26.20 -25.27 -24.55
Mar -20.97 -20.30 -22.38 -21.95 -19.56 -18.90
Apr -15.98 -16.55 -16.14 -16.65 -15.82 -16.55
May -6.89 -7.10 -7.43 -7.30 -6.35 -6.25
Jun 1.66 1.55 1.04 1.00 2.27 2.40
Jul 7.83 7.70 7.27 7.40 8.38 8.40
Aug 7.85 8.05 7.44 8.00 8.26 8.40
Sep 3.57 3.40 3.29 3.05 3.84 3.75
Oct -6.19 -5.65 -6.69 -6.85 -5.68 -4.35
Nov -16.52 -16.40 -16.60 -15.75 -16.44 -16.90
Dec -21.75 -22.10 -21.46 -21.75 -22.03 -22.85
# 9: Create graph
# This graph compares monthly mean air temperature between the two time periods.

ggplot(monthly_summary, aes(x = month, y = mean_temp, fill = period)) +
  geom_col(position = "dodge") +
  scale_x_continuous(
    breaks = 1:12,
    labels = c("Jan", "Feb", "Mar", "Apr", "May", "Jun",
               "Jul", "Aug", "Sep", "Oct", "Nov", "Dec")
  ) +
  labs(
    title = "Monthly Mean Air Temperature by Time Period",
    subtitle = "Seyakha Russian Arctic Weather Station Data",
    x = "Month",
    y = "Mean Air Temperature (°C)",
    fill = "Time Period",
    caption = "Figure 1. This graph compares average monthly air temperature from 1961-1980 and 1981-2000."
  ) +
  theme_minimal()

# 10: Graph with standard deviation
# This graph shows monthly mean temperature with standard deviation bars.

ggplot(monthly_summary, aes(x = month, y = mean_temp, fill = period)) +
  geom_col(position = position_dodge(width = 0.9)) +
  geom_errorbar(
    aes(ymin = mean_temp - sd_temp, ymax = mean_temp + sd_temp),
    position = position_dodge(width = 0.9),
    width = 0.2
  ) +
  scale_x_continuous(
    breaks = 1:12,
    labels = c("Jan", "Feb", "Mar", "Apr", "May", "Jun",
               "Jul", "Aug", "Sep", "Oct", "Nov", "Dec")
  ) +
  labs(
    title = "Monthly Mean Air Temperature with Standard Deviation",
    subtitle = "Seyakha Russian Arctic Weather Station Data",
    x = "Month",
    y = "Mean Air Temperature (°C)",
    fill = "Time Period",
    caption = "Figure 2. This graph shows monthly mean air temperature with standard deviation bars for 1961-1980 and 1981-2000."
  ) +
  theme_minimal()

# 11: Compare monthly mean and median temperatures

mean_median_long <- monthly_summary %>%
  select(period, month, mean_temp, median_temp) %>%
  pivot_longer(
    cols = c(mean_temp, median_temp),
    names_to = "statistic",
    values_to = "temperature"
  ) %>%
  mutate(
    statistic = recode(
      statistic,
      mean_temp = "Mean",
      median_temp = "Median"
    )
  )

ggplot(mean_median_long, aes(x = month, y = temperature, fill = statistic)) +
  geom_col(position = "dodge") +
  facet_wrap(~ period) +
  scale_x_continuous(
    breaks = 1:12,
    labels = c("Jan", "Feb", "Mar", "Apr", "May", "Jun",
               "Jul", "Aug", "Sep", "Oct", "Nov", "Dec")
  ) +
  labs(
    title = "Monthly Mean and Median Air Temperature",
    subtitle = "Seyakha Russian Arctic Weather Station Data",
    x = "Month",
    y = "Air Temperature (°C)",
    fill = "Statistic",
    caption = "Figure 3. This graph compares monthly mean and median air temperature for both time periods."
  ) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

# 12: Statistical test
# This compares the monthly mean temperatures from 1961-1980 and 1981-2000.

period_comparison <- monthly_summary %>%
  select(period, month, mean_temp) %>%
  pivot_wider(
    names_from = period,
    values_from = mean_temp
  ) %>%
  mutate(
    difference = `1981-2000` - `1961-1980`
  )

kable(
  period_comparison,
  caption = "Table 2. Difference in monthly mean temperature between 1961-1980 and 1981-2000."
)
Table 2. Difference in monthly mean temperature between 1961-1980 and 1981-2000.
month 1961-1980 1981-2000 difference
1 -26.710 -24.895 1.815
2 -26.570 -25.270 1.300
3 -22.380 -19.560 2.820
4 -16.140 -15.820 0.320
5 -7.430 -6.355 1.075
6 1.045 2.265 1.220
7 7.270 8.385 1.115
8 7.440 8.255 0.815
9 3.290 3.845 0.555
10 -6.690 -5.685 1.005
11 -16.600 -16.440 0.160
12 -21.460 -22.030 -0.570
temp_t_test <- t.test(
  period_comparison$`1981-2000`,
  period_comparison$`1961-1980`,
  paired = TRUE
)

temp_t_test
## 
##  Paired t-test
## 
## data:  period_comparison$`1981-2000` and period_comparison$`1961-1980`
## t = 3.9398, df = 11, p-value = 0.002313
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
##  0.4277403 1.5105931
## sample estimates:
## mean difference 
##       0.9691667

Conclusion

The results show that the later time period of 1981-2000 generally had warmer monthly mean temperatures than the earlier time period of 1961-1980. The seasonal pattern stayed similar with the coldest temperatures during the winter months and the warmest temperatures during the summer months. However the later period shifted warmer overall.

The paired t-test compares the 12 monthly mean temperatures from both time periods. The p-value was 0.002313, which is below 0.05 so the difference between the two periods is statistically significant. The graphs and summary statistics suggest that the Seyakha Arctic station experienced warming over time.