R Notebook

This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.

Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Ctrl+Shift+Enter.

library(dplyr)  
library(ggplot2)
library(tidyr)  
library(zoo)    

# Load the datasets
df_Tahiti <- read.csv("CE203N_2025 - HW 2 - Data_Tahiti.csv")
df_Darwin <- read.csv("CE203N_2025 - HW 2 - Data_Darwin.csv")

# Ensure the "Year" column is numeric
df_Tahiti$Year <- as.numeric(df_Tahiti$Year)
df_Darwin$Year <- as.numeric(df_Darwin$Year)

# Compute Monthly Pressure Differences (Tahiti - Darwin)
df_Difference <- df_Tahiti
df_Difference[,-1] <- df_Tahiti[,-1] - df_Darwin[,-1]

# Convert to Long Format
df_long <- df_Difference %>%
  pivot_longer(-Year, names_to = "Month", values_to = "Pressure_Diff") 

# Compute Climatology (1951-1980) - Monthly Mean and SD
climatology <- df_long %>%
  filter(Year >= 1951 & Year <= 1980) %>%
  group_by(Month) %>%
  summarise(Mean = mean(Pressure_Diff, na.rm=TRUE), SD = sd(Pressure_Diff, na.rm=TRUE), .groups="drop")

# Compute Standardized SOI
df_long <- df_long %>%
  left_join(climatology, by="Month") %>%
  mutate(SOI = (Pressure_Diff - Mean) / SD) %>%
  select(Year, Month, SOI) %>%
  arrange(Year)

# Convert Month names to numerical order for proper sorting
df_long$Month <- match(df_long$Month, month.abb)  # Convert "Jan" -> 1, "Feb" -> 2, etc.

# Compute 12-Month Running Mean for Smoothing
df_long <- df_long %>%
  arrange(Year, Month) %>%
  mutate(SOI_smoothed = zoo::rollmean(SOI, k = 12, fill = NA, align = "center"))

# Remove NA values
df_long <- df_long %>% filter(!is.na(SOI_smoothed))

# Plot the SOI
ggplot(df_long, aes(x = Year + (Month-1)/12, y = SOI_smoothed, fill = factor(SOI_smoothed > 0))) +
  geom_col(width = 0.5) +  # Reduce width to make bars thinner
  scale_fill_manual(values = c("FALSE" = "red", "TRUE" = "blue"), labels = c("Warm (El Nina)", "Cool (La Nina)")) +
  theme_minimal() +
  scale_x_continuous(breaks = seq(1880, 2020, by = 20)) +  
  labs(title = "Southern Oscillation Index",
       subtitle = "Smoothed 1-Year Running Mean",
       x = "Year", y = "SOI") +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.title = element_blank(),
        axis.text.x = element_text(angle = 45, hjust = 1))

NA
NA

# Constants
T_surf <- 295  # Surface temperature in K
P0 <- 98000  # Surface pressure in Pa
g <- 9.81  # Gravity (m/s^2)
Rd <- 287  # Specific gas constant for dry air (J/kg.K)
Gamma_d <- 9.8 / 1000  # Convert lapse rate to K/m
epsilon <- 0.622  # Ratio of water vapor to dry air

# Generate height range (0 to 5 km)
z <- seq(0, 5000, length.out = 100)  # Height in meters

# Compute temperature profile
T_z <- T_surf - Gamma_d * z

# Compute pressure profile using hydrostatic approximation
P_z <- P0 * (T_z / T_surf)^(g / (Rd * Gamma_d))

# Compute saturation vapor pressure (Clausius-Clapeyron equation)
e_s <- 6.112 * exp((17.67 * (T_z - 273.15)) / (T_z - 29.65)) * 100  # Convert hPa to Pa

# Compute saturated specific humidity
q_s <- (epsilon * e_s) / (P_z - e_s)

# Create a dataframe for ggplot
df <- data.frame(T_z, q_s)

# Plot using ggplot2
ggplot(df, aes(x = z, y = q_s)) +
  geom_line(color = "blue", size = 1.2) +
  labs(x = "Altitude", 
       y = "Saturated Specific Humidity (kg/kg)", 
       title = "Saturated Specific Humidity vs. Altitude") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5, size = 14, face = "bold"))

```