install.packages(‘knitr’)
library(knitr)
print("Welcome to an exploration of Arctic Warming. Here, we will analyze statistics about Agata Station, a taiga ecosystem in Siberia.")[1] "Welcome to an exploration of Arctic Warming. Here, we will analyze statistics about Agata Station, a taiga ecosystem in Siberia."require(ggplot2)Loading required package: ggplot2#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");
# 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
✔ lubridate 1.9.5 ✔ tibble 3.3.1
✔ purrr 1.2.1 ✔ tidyr 1.3.2
── 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[[1]]
[1] TRUE
[[2]]
[1] TRUEinstall.packages("sf") # run once if not installed
The downloaded binary packages are in
/var/folders/dm/sj7012g577qdnv76gj2_szd80000gp/T//RtmpLcizfT/downloaded_packagesinstall.packages("tidyverse")
The downloaded binary packages are in
/var/folders/dm/sj7012g577qdnv76gj2_szd80000gp/T//RtmpLcizfT/downloaded_packageslibrary(sf)Linking to GEOS 3.13.0, GDAL 3.8.5, PROJ 9.5.1; sf_use_s2() is TRUElibrary(tidyverse)# read the file
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")# Compute monthly means
bars <- aggregate(agata[,9], list(t=agata[,3]), mean, na.rm=TRUE)
bars <- rbind(data.frame(t="Ann", x=mean(bars$x)), bars)
# Define ymin
ymin <- floor(min(bars$x, na.rm = TRUE) / 10) * 10
# Shift values so all values are positive
bars_shift <- bars$x - ymin
# Plot
barplot(bars_shift,
yaxt = "n",
names.arg = bars$t,
col = "cadetblue2")
# Create axis starting at 0
ymax <- ceiling(max(bars$x, na.rm = TRUE) / 10) * 10
yticks <- seq(0, ymax - ymin, by = 5)
# Label back in original scale
axis(2, at = yticks, labels = yticks + ymin)
title(main="Mean Temperature",
xlab="Month",
ylab="Temp (°C)")
print("This bar graph shows the monthly mean temperatures at Agata Station for the years 1961 to present.")[1] "This bar graph shows the monthly mean temperatures at Agata Station for the years 1961 to present."df1 <- agata[agata[,2] >= 1961 & agata[,2] <= 1980, ]
df2 <- agata[agata[,2] >= 1981 & agata[,2] <= 2000, ]
bars1 <- aggregate(df1[,9], list(t=df1[,3]), mean, na.rm=TRUE)
bars2 <- aggregate(df2[,9], list(t=df2[,3]), mean, na.rm=TRUE)
bars1 <- rbind(data.frame(t="Ann", x=mean(bars1$x)), bars1)
bars2 <- rbind(data.frame(t="Ann", x=mean(bars2$x)), bars2)
mat <- rbind(bars1$x, bars2$x)
# Find overall minimum across BOTH mean + median
ymin <- floor(min(mat, na.rm = TRUE) / 10) * 10
# Shift data upward
mat_shift <- mat - ymin
bp <- barplot(mat_shift,
beside=TRUE,
names.arg=bars1$t,
col=c("steelblue","cadetblue1"),
legend.text = c("1961–1980", "1981–2000"),
yaxt="n")
yticks <- seq(0, max(mat_shift)+10, by=10)
axis(2, at=yticks, labels=yticks + ymin)
title(main="Temperature Comparison between 1961-1980 and 1981-2000 at Agata Station",
xlab="Month",
ylab="Temp (°C)")
print("This bar graph compares the monthly mean temperatures at Agata Station from the years 1961 to 1980 and the following twenty-year period of and 1981 to 2000.")[1] "This bar graph compares the monthly mean temperatures at Agata Station from the years 1961 to 1980 and the following twenty-year period of and 1981 to 2000."# extract key information; 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 = "cadetblue2",
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 in Agata Station with Standard Deviation",
xlab = "Month",
ylab = "Temperature (Degrees C)")
print("This bar graph compares the monthly mean temperatures at Agata Station from the years 1961 to present, with standard deviation bars present.")[1] "This bar graph compares the monthly mean temperatures at Agata Station from the years 1961 to present, with standard deviation bars present."