# 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."
)
| 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."
)
| 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
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.