knitr::opts_chunk$set(echo = TRUE)Agata and Indiga: A Comparative Study
The Changing Climate of Siberia: An Analysis of The Impacts of Anthropogenic Climate Change on Northern Russia
#Loading the Packages
options(repos = c(CRAN = "https://cloud.r-project.org"))
# Load a list of packages. Install them first if they are not available.
# The list of packages to be installed
list.of.packages <- c("tidyverse", "ggplot2","maps","ggthemes");
# Check out the packages that have not been installed yet.
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
# Install those missing packages first. It could take a long time for the first time.
if(length(new.packages)>0) install.packages(new.packages)
# Load all packages.
lapply(list.of.packages,function(x) {
require(x,character.only = TRUE,quietly = TRUE)
})── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr 1.2.0 ✔ readr 2.1.6
✔ forcats 1.0.1 ✔ stringr 1.6.0
✔ ggplot2 4.0.3 ✔ tibble 3.3.1
✔ lubridate 1.9.5 ✔ tidyr 1.3.2
✔ purrr 1.2.1
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag() masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
Attaching package: 'maps'
The following object is masked from 'package:purrr':
map[[1]]
[1] TRUE
[[2]]
[1] TRUE
[[3]]
[1] TRUE
[[4]]
[1] TRUEMap of the Locations of Indiga and Agata Stations
install.packages("ggplot2")
The downloaded binary packages are in
/var/folders/dm/sj7012g577qdnv76gj2_szd80000gp/T//RtmpY1FBkM/downloaded_packagesinstall.packages("maps")
The downloaded binary packages are in
/var/folders/dm/sj7012g577qdnv76gj2_szd80000gp/T//RtmpY1FBkM/downloaded_packagesinstall.packages("ggthemes")
The downloaded binary packages are in
/var/folders/dm/sj7012g577qdnv76gj2_szd80000gp/T//RtmpY1FBkM/downloaded_packageslibrary(ggthemes)
library(ggplot2)
library(maps)
# Coordinates sourced from the data
indigadf <- data.frame(
station = "Indiga",
lon = 93.47,
lat = 66.88
)
agatadf <- data.frame(
station = "Agata",
lon = 48.77,
lat = 66.70
)
# Combine into one dataframe
df <- rbind(indigadf, agatadf)
# World map
world <- map_data("world")
# Plot
ggplot() +
geom_polygon(data = world,
aes(x = long, y = lat, group = group),
fill = "darkseagreen",
color = "white") +
geom_point(data = df,
aes(x = lon, y = lat),
color = "blue",
size = 4) +
geom_text(data = df,
aes(x = lon, y = lat, label = station),
hjust = -0.2,
vjust = 0.5) +
coord_fixed(1.3,
xlim = c(20, 150),
ylim = c(40, 80)) +
labs(title = "Locations of Indiga and Agata Stations",
x = "Longitude",
y = "Latitude") +
theme_clean() +
theme(
plot.title = element_text(hjust = 0.5)
)
Temperature Comparision
url <- "https://noaadata.apps.nsidc.org/NOAA/G02141/uni.agata.23383.dat"
download.file(url,destfile='data.dat',method='curl')
# Read text file (999.99 is NA)
agata <- read.table('data.dat',sep="",na.strings="999.99")
url <- "https://noaadata.apps.nsidc.org/NOAA/G02141/uni.indiga.22292.dat"
download.file(url,destfile='data.dat',method='curl')
# Read text file (999.99 is NA)
indiga <- read.table('data.dat',sep="",na.strings="999.99")
# Make sure to set up a matrix first
layout(matrix(c(1,2), nrow = 1, ncol = 2))
# create a graphing function
plot_temp <- function(dat, station_name){
# Compute monthly means
bars <- aggregate(dat[,9], list(t = dat[,3]), mean, na.rm = TRUE)
bars <- rbind(data.frame(t = "Ann", x = mean(bars$x)), bars)
# Define ymax and ymin limits
ymin <- floor(min(bars$x, na.rm = TRUE) / 10) * 10
ymax <- ceiling(max(bars$x, na.rm = TRUE) / 10) * 10
# Shift values so minimum becomes zero
bars_shift <- bars$x - ymin
# Tick marks
yticks <- seq(0, ymax - ymin, by = 5)
# Plot
barplot(
bars_shift,
yaxt = "n",
ylim = c(0, ymax - ymin),
names.arg = bars$t,
col = "cadetblue2",
border = "black"
)
# Custom y-axis
axis(2, at = yticks, labels = yticks + ymin)
# Titles
title(
main = paste("Mean Temperature at", station_name),
xlab = "Month",
ylab = "Temp (°C)"
)
}
# Plot both stations side by side
plot_temp(agata, "Agata Station")
plot_temp(indiga, "Indiga Station")
Temperature Over Time
url <- "https://noaadata.apps.nsidc.org/NOAA/G02141/uni.agata.23383.dat"
download.file(url,destfile='data.dat',method='curl')
# Read text file (999.99 is NA for Temperature)
agata <- read.table('data.dat',sep="",na.strings="999.99")
url <- "https://noaadata.apps.nsidc.org/NOAA/G02141/uni.indiga.22292.dat"
download.file(url,destfile='data.dat',method='curl')
# Read text file (999.99 is NA for Temperature)
indiga <- read.table('data.dat',sep="",na.strings="999.99")
#Temperature Data from Indiga
# create data frames x 2
df1 <- indiga[indiga[,2] >= 1961 & indiga[,2] <= 1980, ]
df2 <- indiga[indiga[,2] >= 1981 & indiga[,2] <= 2000, ]
# make bars
bars1 <- aggregate(df1[,9], list(t=df1[,3]), mean, na.rm=TRUE)
bars2 <- aggregate(df2[,9], list(t=df2[,3]), mean, na.rm=TRUE)
# add annual bar
bars1 <- rbind(data.frame(t="Ann", x=mean(bars1$x)), bars1)
bars2 <- rbind(data.frame(t="Ann", x=mean(bars2$x)), bars2)
# Monthly Means Between 1961-1980 and 1981-2000: Agata Station
# create data frames
df3 <- agata[agata[,2] >= 1961 & agata[,2] <= 1980, ]
df4 <- agata[agata[,2] >= 1981 & agata[,2] <= 2000, ]
# make bars
bars3 <- aggregate(df3[,19], list(t=df3[,3]), mean, na.rm=TRUE)
bars4 <- aggregate(df4[,19], list(t=df4[,3]), mean, na.rm=TRUE)
# add annual bar
bars3 <- rbind(data.frame(t="Ann", x=mean(bars3$x)), bars3)
bars4 <- rbind(data.frame(t="Ann", x=mean(bars4$x)), bars4)
# matrix for grouped bars
mat <- rbind(bars3$x, bars4$x)
#Creating a matrix to display side-by-side comparisons
layout(matrix(c(1,2), nrow = 1, ncol = 2))
# matrix for grouped bars
mat <- rbind(bars1$x, bars2$x)
ymin <- -40
mat_shift <- mat - ymin
ymax <- ceiling(max(mat, na.rm = TRUE)/10)*10
yticks <- seq(0, ymax - ymin, by = 10)
barplot(mat_shift,
beside = TRUE,
names.arg = bars1$t,
ylim = c(0, ymax - ymin),
yaxt = "n",
col = c("cyan4", "cyan3"),
legend.text = c("1961–1980", "1981–2000"),
args.legend = list(x = "topleft"))
axis(2, at = yticks, labels = yticks + ymin)
title(main="Temperature Comparison: Indiga Station",
xlab="Month",
ylab="Temperature in °C")
# matrix for grouped bars
mat <- rbind(bars3$x, bars4$x)
ymin <- -40
mat_shift <- mat - ymin
ymax <- ceiling(max(mat, na.rm = TRUE)/10)*10
yticks <- seq(0, ymax - ymin, by = 10)
barplot(mat_shift,
beside = TRUE,
names.arg = bars3$t,
ylim = c(0, ymax - ymin),
yaxt = "n",
col = c("cyan4", "cyan3"),
legend.text = c("1961–1980", "1981–2000"),
args.legend = list(x = "topleft"))
axis(2, at = yticks, labels = yticks + ymin)
title(main="Temperature Comparison: Agata Station",
xlab="Month",
ylab="Temperature in °C")
Standard Deviation
# extract key information about Agata Station; find the standard deviation
bars <- aggregate(agata[,9], list(t=agata[,3]), mean, na.rm=TRUE)
sds <- aggregate(agata[,9], list(t=agata[,3]), sd, na.rm=TRUE)
# do the same to add an annual bar with standard deviation
bars <- rbind(data.frame(t="Ann", x=mean(bars$x)), bars)
sds <- rbind(data.frame(t="Ann", x=mean(sds$x)), sds)
# shift
shift <- min(bars$x - sds$x)
bars$x <- bars$x - shift
sds$x <- sds$x
# calculate ymax
ymax <- max(bars$x + sds$x)
# Round up to nearest 5
ymax <- ceiling(ymax / 5) * 5
ymax <- round(ymax)
# Plot
bp <- barplot(bars$x,
names.arg = bars$t,
col = "deepskyblue3",
yaxt = "n",
ylim = c(0, ymax))
# Add the error bars
arrows(bp,
bars$x - sds$x,
bp,
bars$x + sds$x,
angle = 90, code = 3, length = 0.05)
# Make an axis
yticks <- seq(0, ymax, by = 5)
axis(2,
at = yticks,
labels = round(yticks + shift, 0))
# Add the title and labels
title(main = "Temperature at Agata Station with Standard Deviation",
xlab = "Month",
ylab = "Temperature in °C")
# extract key information about Indiga Staion; find the standard deviation
bars1 <- aggregate(indiga[,9], list(t=indiga[,3]), mean, na.rm=TRUE)
sds1 <- aggregate(indiga[,9], list(t=indiga[,3]), sd, na.rm=TRUE)
# do the same to add an annual bar with standard deviation
bars1 <- rbind(data.frame(t="Ann", x=mean(bars1$x)), bars1)
sds1 <- rbind(data.frame(t="Ann", x=mean(sds1$x)), sds1)
# shift
shift <- min(bars1$x - sds1$x)
bars1$x <- bars1$x - shift
sds1$x <- sds1$x
# calculate ymax
ymax <- max(bars1$x + sds1$x)
# Round up to nearest 5
ymax <- ceiling(ymax / 5) * 5 +5
ymax <- round(ymax)
# Plot
bp <- barplot(bars1$x,
names.arg = bars1$t,
col = "deepskyblue3",
yaxt = "n",
ylim = c(0, ymax))
# Add the error bars
arrows(bp,
bars1$x - sds1$x,
bp,
bars1$x + sds1$x,
angle = 90, code = 3, length = 0.05)
# Make an axis
yticks <- seq(0, ymax, by = 5)
axis(2,
at = yticks,
labels = round(yticks + shift, 0))
# Add the title and labels
title(main = "Temperature at Indiga Station with Standard Deviation",
xlab = "Month",
ylab = "Temperature in °C")
Statistical Analysis
#plot year against temperature at Indiga Station
plot(indiga$V2, indiga$V9,
pch = 16,
col = rgb(0, 0.1, 0.2, 0.3),
xlab = "Year",
ylab = "Temperature",
main = "Indiga Temperatures with Line of Best Fit")
# create linear model
model <- lm(indiga$V9 ~ indiga$V2)
# create a line of best fit
abline(model, col = "red", lwd = 2)
#plot year against temperature at Agata Station
plot(agata$V2, agata$V9,
pch = 16,
col = rgb(0, 0.1, 0.2, 0.3),
xlab = "Year",
ylab = "Temperature",
main = "Agata Temperatures with Line of Best Fit")
# create a line of best fit
# create linear model
model <- lm(agata$V9 ~ agata$V2)
abline(model, col = "red", lwd = 2)
Precipation Plots
#Download Agata Station Data
url <- "https://noaadata.apps.nsidc.org/NOAA/G02141/uni.agata.23383.dat"
download.file(url,destfile='data.dat',method='curl')
# Read text file (-1.00 is NA for precip)
agata <- read.table('data.dat',sep="",na.strings="-1.00")
#Download Indiga Station Data
url <- "https://noaadata.apps.nsidc.org/NOAA/G02141/uni.indiga.22292.dat"
download.file(url,destfile='data.dat',method='curl')
# Read text file (-1.00 is NA for precip)
indiga <- read.table('data.dat',sep="",na.strings="-1.00")
layout(matrix(c(1,2), nrow = 1, ncol = 2))
# plot
barplot(mat,
beside = TRUE,
names.arg = bars3$t,
ylim = c(0, 70),
col = c("cyan4", "cyan3"),
legend.text = c("1961–1980", "1981–2000"),
args.legend = list(x = "topleft"))
title(main="Precipitation Comparison: Agata Station",
xlab="Month",
ylab="Precipitation in Millimeters")
# Monthly Means Between 1961-1980 and 1981-2000 : Precipitation at Indiga Station
# create data frames x 2
df1 <- indiga[indiga[,2] >= 1961 & indiga[,2] <= 1980, ]
df2 <- indiga[indiga[,2] >= 1981 & indiga[,2] <= 2000, ]
# make bars
bars1 <- aggregate(df1[,19], list(t=df1[,3]), mean, na.rm=TRUE)
bars2 <- aggregate(df2[,19], list(t=df2[,3]), mean, na.rm=TRUE)
# add annual bar
bars1 <- rbind(data.frame(t="Ann", x=mean(bars1$x)), bars1)
bars2 <- rbind(data.frame(t="Ann", x=mean(bars2$x)), bars2)
# make a matrix to bind the two bar graphs
mat <- rbind(bars1$x, bars2$x)
barplot(mat,
beside = TRUE,
names.arg = bars1$t,
ylim = c(0, 70),
col = c("cyan4", "cyan3"),
legend.text = c("1961–1980", "1981–2000"),
args.legend = list(x = "topleft"))
title(main="Precipitation Comparison: Indiga Station",
xlab="Month",
ylab="Precipitation in Millimeters")