Post-Homogenization Temperature Analysis

Introduction

This tutorial guides you through the post-homogenization analysis of temperature data, including visualization and anomaly calculation.

Step 1: Load Required Libraries

We begin by loading the necessary libraries for data processing and visualization.

# Load required libraries
packages <- c("dplyr", "tidyr", "lubridate", "ggplot2", "dygraphs", "xts")

# Install missing packages
installed_packages <- packages %in% rownames(installed.packages())
if (any(installed_packages == FALSE)) {
  install.packages(packages[!installed_packages])
}

# Load the packages
invisible(lapply(packages, library, character.only = TRUE))

Step 2: Read and Transform Homogenized Data

We read the output from the homogenization process and convert it into a long format.

tn.hom <- read.csv("tn_1981-2022_series.csv")

tn.hom <- tn.hom |> 
         pivot_longer(-Date, values_to = "TN_hom", names_to = "SOUID") |> 
         mutate(SOUID = as.character(gsub("X", "", SOUID)),
                Date = ymd(Date))

Step 3: Filter and Aggregate Data

### Filter one station
tn.hom.f <- tn.hom |> filter(SOUID == "100682")

### Calculate monthly means
tn.hom.fm <- tn.hom |> 
              group_by(format(Date, "%Y-%m")) |> 
              summarise(Mn = mean(TN_hom)) |> 
              rename("Date" = `format(Date, "%Y-%m")`)

Step 4: Interactive Time Series Visualization

We use dygraphs to visualize the temperature trends interactively.

dy <- xts(x = tn.hom.f$TN_hom, order.by = tn.hom.f$Date)

p <- dygraph(dy) %>%
  dyOptions(labelsUTC = TRUE, fillGraph=TRUE, fillAlpha=0.1, drawGrid = FALSE, colors="darkgrey") %>%
  dyRangeSelector() %>%
  dyCrosshair(direction = "vertical") %>%
  dyHighlight(highlightCircleSize = 5, highlightSeriesBackgroundAlpha = 0.2, hideOnMouseOut = FALSE) %>%
  dyRoller(rollPeriod = 1)

p

Step 5: Anomaly Calculations

### Calculate monthly anomalies
tn.monmean <- tn.hom |> group_by(SOUID, format(Date, "%Y-%m")) |> 
  summarise(Tn_mn = mean(TN_hom)) |> 
  mutate(Date = as.Date(paste0(`format(Date, "%Y-%m")`, "-01"))) |> 
  group_by(SOUID, month(Date)) |> 
  mutate(Tn_ymn = mean(Tn_mn),
         Anom = Tn_mn - Tn_ymn,
         Year = year(Date),
         Month = month(Date))

Step 6: Compute Station and Overall Temperature Anomalies

## make a mean across all the station 
tn.monavg <- tn.hom |> group_by(format(Date,"%Y-%m")) |> 
  summarise(Tn_mn = mean(TN_hom)) |> 
  mutate(Date = as.Date(paste0(`format(Date, "%Y-%m")`,"-01"))) |> 
  group_by(month(Date)) |> 
  mutate(Tn_ymn = mean(Tn_mn),
         Anom = Tn_mn - Tn_ymn,
         Year = year(Date),
         Month = month(Date))

Step 7: Yearly Anomalies

tn.ymean <- tn.hom |> group_by(SOUID, format(Date, "%Y")) |> 
  summarise(Tn_mn = mean(TN_hom)) |> 
  group_by(SOUID) |> 
  mutate(Tn_ymn = mean(Tn_mn),
         Anom = Tn_mn - Tn_ymn) |> 
  rename("Date"= `format(Date, "%Y")`)
tn.ymean$Date <- as.integer(tn.ymean$Date)

## make a mean across all the station 
tn.yearavg <- tn.hom|>group_by(format(Date,"%Y"))|>
  summarise(Tn_mn = mean(TN_hom))|>
  mutate(Tn_ymn = mean(Tn_mn),
         Anom = Tn_mn - Tn_ymn)|>
  rename( "Date"= `format(Date, "%Y")`)

tn.yearavg$Date <- as.integer(tn.yearavg$Date)

Step 8: Visualization of Anomalies

Yearly Anomaly Plot

ggplot(tn.yearavg) +
  geom_col(aes(x = Date, y = Anom, fill = Anom), show.legend = FALSE) +
  scale_fill_gradient2(low="#62b2e4", high="#ab0e14", mid="grey", midpoint=0) +
  scale_x_continuous(breaks = seq(1981, 2022, 5)) +
  theme(axis.text.x = element_text(angle = 70, vjust = .1, hjust = .3)) +
  xlab("") + ylab("Anomaly [°C]") + theme_bw()

Yearly Absolute Values Plot

#### lines for absolute values
ggplot(tn.yearavg,aes(
  x = Date,
  y = Tn_mn)) +
  geom_line(size=.4, show.legend = F) +
  geom_point(color = "black", show.legend = F)+
  geom_smooth(color = "red", se = F, method = "lm")+
  scale_x_continuous(breaks =  seq(1981,2022,5))+
  theme(axis.text.x =  element_text(angle = 70, vjust = .1, hjust = .3))+
  xlab("")+ylab(" Air temperature [°C]")+theme_bw()

Montlhy Anomaly Plot

#### lines for absolute values
ggplot(tn.monavg) +
  geom_col(aes(
    x = Date,
    y = Anom,
    fill = Anom
  ), show.legend = F) +
  scale_fill_gradient2(low="#62b2e4", high="#ab0e14", mid="grey",midpoint= 0)+
  scale_x_continuous(breaks =  seq(1981,2022,5))+
  theme(axis.text.x =  element_text(angle = 70, vjust = .1, hjust = .3))+
  xlab("")+ylab("Anomaly [°C]")+theme_bw()+facet_wrap(~Month)

Monthly Absolute Plot

#### lines for absolute values
ggplot(tn.monavg,aes(
  x = Date,
  y = Tn_mn)) +
  geom_line(size=.4, show.legend = F) +
  geom_point(color = "black", show.legend = F)+
  geom_smooth(color = "red", se = F, method = "lm")+
  scale_x_continuous(breaks =  seq(1981,2022,5))+
  theme(axis.text.x =  element_text(angle = 70, vjust = .1, hjust = .3))+
  xlab("")+ylab(" Air temperature [°C]")+theme_bw()+facet_wrap(~Month, scales = "free")

Step 9:Station-Level Visualization

### let's plot each station one row 
ggplot(data = tn.ymean)+
        geom_tile(aes(x = Date, y = SOUID, fill =  Anom))+
      scale_fill_gradient2(low="#62b2e4", high="#ab0e14", mid="white",midpoint= 0)+
  scale_x_continuous(breaks =  seq(1981,2022,5), expand = c(0,0))+
  theme(axis.text.x =  element_text(angle = 70, vjust = .1, hjust = .3))+
  xlab("")+ylab("Station's ID")+theme_bw()

Step 10: Monthly Station-Level Visualization

 ggplot(data = tn.monmean)+
  geom_raster(aes(x = Date, y = SOUID, fill =  Anom))+
  scale_fill_gradient2(low="#62b2e4", high="#ab0e14", mid="white",midpoint= 0)+
  scale_x_continuous(breaks =  seq(1981,2022,5), expand = c(0,0))+
  theme(axis.text.x =  element_text(angle = 70, vjust = .1, hjust = .3))+
  xlab("")+ylab("Station's ID")+theme_bw()+facet_wrap(~Month)

Conclusion

This tutorial covers the post-homogenization processing of temperature data, including anomaly calculations and visualizations.