Heatwave/Temperature Analysis Using API call from nasapower package

Author

Ukpai Chinedu Kalu

API and LIBRARY NEEDED

# Install if not yet installed
library(nasapower)
library(dplyr)
library(flextable)
library(lubridate)
library(ggplot2)
library(patchwork)

DATA PREP

# Coordinates for Nsukka, Enugu State
lat <- 6.8562
lon <- 7.3958

# Download daily temperature data (1981 - present)
temp_data <- get_power(
  community = "AG",
  lonlat = c(lon, lat),
  pars = c("T2M_MAX", "T2M_MIN", "T2M"), # max, min, avg temp
  dates = c("1981-01-01", "2023-12-31"),
  temporal_api = "DAILY"
)

# View data
head(temp_data)

────────────────────────────────────────────────────────────────────────────────

── NASA/POWER Source Native Resolution Daily Data  ─────────────────────────────

Dates (month/day/year): 01/01/1981 through 12/31/2023 in LST

Location: latitude 6.8562 longitude 7.3958

elevation from MERRA-2: Average for 0.5 x 0.625 degree lat/lon region = 220.37
meters

The value for missing source data that cannot be computed or is outside of the
sources availability range: NA

parameter(s):

────────────────────────────────────────────────────────────────────────────────

Parameters:

T2M_MAX MERRA-2 Temperature at 2 Meters Maximum (C) ; T2M_MIN MERRA-2
Temperature at 2 Meters Minimum (C) ; T2M MERRA-2 Temperature at 2 Meters (C)

────────────────────────────────────────────────────────────────────────────────
# A tibble: 6 × 10
    LON   LAT  YEAR    MM    DD   DOY YYYYMMDD   T2M_MAX T2M_MIN   T2M
  <dbl> <dbl> <dbl> <int> <int> <int> <date>       <dbl>   <dbl> <dbl>
1  7.40  6.86  1981     1     1     1 1981-01-01    30.3    18.7  24.5
2  7.40  6.86  1981     1     2     2 1981-01-02    29.8    18.3  23.7
3  7.40  6.86  1981     1     3     3 1981-01-03    27.8    17.2  22.0
4  7.40  6.86  1981     1     4     4 1981-01-04    26.2    16.0  20.6
5  7.40  6.86  1981     1     5     5 1981-01-05    27.7    14.7  20.2
6  7.40  6.86  1981     1     6     6 1981-01-06    27      14.3  20.2

# Convert to data frame
temp_df <- as.data.frame(temp_data)

# Quick plot
plot(temp_df$T2M ~ as.Date(temp_df$YYYYMMDD),
     type="l", col="red",
     main="Daily Average Temperature in Nsukka",
     xlab="Date", ylab="Temperature (°C)")

Clean the Data Frame

The raw NASA POWER output has many metadata columns you may not need.

keep only the date + temperature variables:

  Year DayOfYear       Date Temp_Avg Temp_Min Temp_Max
1 1981         1 1981-01-01    24.48    18.66    30.31
2 1981         2 1981-01-02    23.67    18.31    29.79
3 1981         3 1981-01-03    21.99    17.20    27.83
4 1981         4 1981-01-04    20.56    15.96    26.23
5 1981         5 1981-01-05    20.19    14.73    27.70
6 1981         6 1981-01-06    20.23    14.32    27.00

Export to CSV

write.csv(temp_pub, "Nsukka_Temperature.csv", row.names = FALSE)

Export a Formatted Table (for direct use in papers)

# Take a sample (first 10 rows) for publication table
temp_sample <- head(temp_pub, 10)

# Create a nice table
flextable(temp_sample)

Year	DayOfYear	Date	Temp_Avg	Temp_Min	Temp_Max
1,981	1	1981-01-01	24.48	18.66	30.31
1,981	2	1981-01-02	23.67	18.31	29.79
1,981	3	1981-01-03	21.99	17.20	27.83
1,981	4	1981-01-04	20.56	15.96	26.23
1,981	5	1981-01-05	20.19	14.73	27.70
1,981	6	1981-01-06	20.23	14.32	27.00
1,981	7	1981-01-07	20.93	16.27	26.74
1,981	8	1981-01-08	22.58	15.90	28.81
1,981	9	1981-01-09	23.80	19.81	28.59
1,981	10	1981-01-10	22.63	19.19	27.74

Summarize Monthly Data

temp_pub <- temp_pub %>%
  mutate(Year = year(Date),
         Month = month(Date, label = TRUE, abbr = TRUE))

# Monthly summary (average, min, max)
monthly_summary <- temp_pub %>%
  group_by(Year, Month) %>%
  summarise(
    Avg_Temp = mean(Temp_Avg, na.rm = TRUE),
    Min_Temp = mean(Temp_Min, na.rm = TRUE),
    Max_Temp = mean(Temp_Max, na.rm = TRUE),
    .groups = "drop"
  )

head(monthly_summary)

# A tibble: 6 × 5
   Year Month Avg_Temp Min_Temp Max_Temp
  <dbl> <ord>    <dbl>    <dbl>    <dbl>
1  1981 Jan       23.0     17.8     28.7
2  1981 Feb       24.8     19.2     30.6
3  1981 Mar       25.7     22.1     30.3
4  1981 Apr       26.2     22.9     30.4
5  1981 May       25.1     22.4     28.7
6  1981 Jun       25.1     22.3     28.6

Summarize Annual Data

annual_summary <- temp_pub %>%
  group_by(Year) %>%
  summarise(
    Avg_Temp = mean(Temp_Avg, na.rm = TRUE),
    Min_Temp = mean(Temp_Min, na.rm = TRUE),
    Max_Temp = mean(Temp_Max, na.rm = TRUE),
    .groups = "drop"
  )

annual_summary

# A tibble: 43 × 4
    Year Avg_Temp Min_Temp Max_Temp
   <dbl>    <dbl>    <dbl>    <dbl>
 1  1981     24.4     20.8     28.7
 2  1982     24.1     20.7     28.2
 3  1983     24.1     20.4     28.6
 4  1984     24.1     20.5     28.3
 5  1985     24.3     20.6     28.6
 6  1986     24.4     20.7     28.9
 7  1987     25.3     21.5     29.9
 8  1988     25.0     21.5     29.2
 9  1989     24.0     20.3     28.6
10  1990     25.0     21.3     29.4
# ℹ 33 more rows

# Table for publication
flextable(annual_summary)

Year	Avg_Temp	Min_Temp	Max_Temp
1,981	24.38088	20.78995	28.68885
1,982	24.10444	20.66164	28.22151
1,983	24.12964	20.42222	28.56992
1,984	24.06683	20.46921	28.33046
1,985	24.27277	20.61077	28.57008
1,986	24.38956	20.74386	28.90723
1,987	25.34729	21.50449	29.93518
1,988	24.99992	21.47877	29.24385
1,989	24.03600	20.26211	28.58995
1,990	25.01753	21.34027	29.37460
1,991	25.02838	21.38556	29.65074
1,992	24.69552	20.73801	29.60566
1,993	24.49304	20.93671	28.80874
1,994	24.34532	20.85721	28.59559
1,995	24.49345	21.07156	28.60507
1,996	24.63801	21.13672	28.91022
1,997	24.55573	20.93655	28.88792
1,998	25.18290	21.45329	29.63227
1,999	25.35512	21.60786	30.07625
2,000	25.22246	21.15281	30.22063
2,001	24.52101	20.74058	28.93636
2,002	24.72268	21.02151	29.13978
2,003	25.07090	21.54175	29.35066
2,004	25.04683	21.57754	29.40202
2,005	25.49663	21.90373	30.01463
2,006	25.37066	21.82151	29.75307
2,007	24.71611	21.16504	29.05458
2,008	24.70322	21.06533	28.91702
2,009	25.21293	21.77290	29.31405
2,010	25.26841	21.80014	29.40274
2,011	24.66312	21.08605	28.81085
2,012	24.74197	21.36866	28.81112
2,013	24.99301	21.69751	29.05912
2,014	25.31403	21.80195	29.67184
2,015	25.60704	21.77477	30.48912
2,016	26.00563	22.10738	30.79948
2,017	25.75044	22.08981	30.23030
2,018	25.43512	21.71060	29.88904
2,019	25.70079	22.01123	30.05416
2,020	25.52481	21.57145	30.23128
2,021	25.65934	21.89830	30.30633
2,022	25.16452	21.24192	30.08556
2,023	26.14282	22.21178	30.97605

Export Tables for Publication

# Export monthly summary
write.csv(monthly_summary, "Nsukka_Temperature_Monthly.csv", row.names = FALSE)

# Export annual summary
write.csv(annual_summary, "Nsukka_Temperature_Annual.csv", row.names = FALSE)

Define Extreme Temperature

For heatwave detection, we usually need a threshold.
Common definitions are:

Percentile-based: days above the 90th or 95th percentile of Tmax.
Absolute threshold: e.g., Tmax > 35 °C for several days.
We are using the 95th percentile of Tmax for this study.

# 1. Define heatwave threshold (95th percentile of Tmax)
threshold <- quantile(temp_pub$Temp_Max, 0.95, na.rm = TRUE)

# Flag extreme days
temp_pub <- temp_pub %>%
  mutate(Extreme = Temp_Max >= threshold,
         Year = year(Date),
         Month = month(Date, label = TRUE, abbr = TRUE))

# --- 2. Monthly Summary ---
heatwave_monthly <- temp_pub %>%
  filter(Extreme) %>%
  group_by(Year, Month) %>%
  summarise(Extreme_Days = n(), .groups = "drop")

# --- 3. Heatmap of Extreme-Day Counts ---
ggplot(heatwave_monthly, aes(x = Year, y = Month, fill = Extreme_Days)) +
  geom_tile(color = "white") +
  scale_fill_gradient(low = "yellow", high = "red") +
  labs(title = "Heatwave Occurrence in Nsukka",
       x = "Year", y = "Month",
       fill = "Extreme Days") +
  theme_minimal()

# --- 4. Consecutive Heatwave Events (>= 3 days) ---
temp_pub <- temp_pub %>%
  arrange(Date) %>%
  mutate(Heatwave_Group = cumsum(!Extreme))  # group consecutive days

heatwave_events <- temp_pub %>%
  filter(Extreme) %>%
  group_by(Heatwave_Group) %>%
  summarise(
    Start = min(Date),
    End = max(Date),
    Duration = n(),
    Max_Temp = max(Temp_Max, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  filter(Duration >= 3)  # only keep 3+ day heatwaves

heatwave_events

# A tibble: 83 × 5
   Heatwave_Group Start      End        Duration Max_Temp
            <int> <date>     <date>        <int>    <dbl>
 1            803 1983-03-17 1983-03-20        4     34.2
 2           2233 1987-03-02 1987-03-06        5     34.6
 3           2260 1987-04-07 1987-04-11        5     34.4
 4           2269 1987-04-22 1987-04-24        3     34.9
 5           2270 1987-04-26 1987-04-30        5     34.6
 6           2562 1988-02-19 1988-02-22        4     34.2
 7           2564 1988-02-25 1988-02-28        4     34.5
 8           3618 1991-01-29 1991-02-06        9     33.9
 9           3623 1991-02-14 1991-02-16        3     34.1
10           3952 1992-01-23 1992-02-03       12     35.5
# ℹ 73 more rows

# --- 5. Max Temperature by Month ---
max_months <- temp_pub %>%
  group_by(Year, Month) %>%
  summarise(Max_Temp = max(Temp_Max, na.rm = TRUE), .groups = "drop")

ggplot(max_months, aes(x = Year, y = Month, fill = Max_Temp)) +
  geom_tile(color = "white") +
  scale_fill_gradient(low = "lightblue", high = "darkred") +
  labs(title = "Extreme Maximum Monthly Temperatures - Nsukka",
       x = "Year", y = "Month",
       fill = "Max Temp (°C)") +
  theme_minimal()

Heatwave Map (Calendar Heatmap Style)

Heatmap of months vs years showing frequency of extreme days.

ggplot(heatwave_monthly, aes(x = Year, y = Month, fill = Extreme_Days)) +
  geom_tile(color = "white") +
  scale_fill_gradient(low = "yellow", high = "red") +
  labs(title = "Heatwave Occurrence in Nsukka",
       x = "Year", y = "Month",
       fill = "Extreme Days") +
  theme_minimal()

Highlight Maximum Heatwave Months

Months with the absolute highest Tmax values;

max_months <- temp_pub %>%
  group_by(Year, Month = month(Date, label = TRUE, abbr = TRUE)) %>%
  summarise(Max_Temp = max(Temp_Max, na.rm = TRUE), .groups = "drop")

ggplot(max_months, aes(x = Year, y = Month, fill = Max_Temp)) +
  geom_tile(color = "white") +
  scale_fill_gradient(low = "lightblue", high = "darkred") +
  labs(title = "Months of Extreme Maximum Temperatures (Nsukka)",
       x = "Year", y = "Month",
       fill = "Max Temp (°C)") +
  theme_minimal()

Export to CSV

write.csv(heatwave_events, "Nsukka_Heatwave_Events.csv", row.names = FALSE)

Create Publication-Ready Table

Formatted table for your paper (APA/journal style):

flextable(heatwave_events) %>%
  set_caption("Table: Detected Heatwave Events in Nsukka") %>%
  autofit()

Heatwave_Group	Start	End	Duration	Max_Temp
803	1983-03-17	1983-03-20	4	34.24
2,233	1987-03-02	1987-03-06	5	34.60
2,260	1987-04-07	1987-04-11	5	34.45
2,269	1987-04-22	1987-04-24	3	34.87
2,270	1987-04-26	1987-04-30	5	34.62
2,562	1988-02-19	1988-02-22	4	34.17
2,564	1988-02-25	1988-02-28	4	34.49
3,618	1991-01-29	1991-02-06	9	33.88
3,623	1991-02-14	1991-02-16	3	34.06
3,952	1992-01-23	1992-02-03	12	35.53
3,953	1992-02-05	1992-02-23	19	36.97
3,955	1992-02-26	1992-03-08	12	38.16
3,957	1992-03-13	1992-03-15	3	35.13
6,103	1998-02-04	1998-02-08	5	34.20
6,133	1998-03-15	1998-03-17	3	33.89
6,137	1998-03-23	1998-03-27	5	34.87
6,138	1998-03-29	1998-03-31	3	35.97
6,491	1999-03-24	1999-03-26	3	34.25
6,750	1999-12-18	1999-12-23	6	34.22
6,756	1999-12-31	2000-01-02	3	34.33
6,762	2000-01-11	2000-01-18	8	36.53
6,763	2000-01-20	2000-01-24	5	35.53
6,764	2000-01-26	2000-01-31	6	36.44
6,766	2000-02-03	2000-02-05	3	34.30
6,769	2000-02-11	2000-03-25	44	37.57
6,771	2000-03-28	2000-04-01	5	37.23
6,784	2000-04-25	2000-05-08	14	36.42
6,785	2000-05-10	2000-05-13	4	37.19
7,438	2002-03-03	2002-03-06	4	34.19
8,163	2004-03-13	2004-03-15	3	33.66
8,165	2004-03-19	2004-03-22	4	34.17
8,478	2005-02-01	2005-02-04	4	34.39
8,821	2006-01-19	2006-01-22	4	33.77
8,833	2006-02-06	2006-02-10	5	35.15
8,886	2006-04-10	2006-04-12	3	34.88
9,220	2007-03-17	2007-03-20	4	34.33
10,312	2010-03-28	2010-03-30	3	34.23
12,063	2015-01-20	2015-01-24	5	35.46
12,067	2015-01-30	2015-02-05	7	35.66
12,084	2015-02-28	2015-03-08	9	35.29
12,087	2015-03-12	2015-03-15	4	34.25
12,115	2015-04-14	2015-04-19	6	34.91
12,332	2015-11-28	2015-11-30	3	33.91
12,358	2015-12-27	2015-12-30	4	33.88
12,361	2016-01-03	2016-01-24	22	36.09
12,365	2016-01-29	2016-03-07	39	37.87
12,368	2016-03-14	2016-03-17	4	35.01
12,376	2016-03-28	2016-03-30	3	34.87
12,388	2016-04-14	2016-04-17	4	34.48
12,655	2017-01-16	2017-01-22	7	34.80
12,663	2017-02-03	2017-02-07	5	34.68
12,664	2017-02-09	2017-02-14	6	35.24
12,669	2017-02-20	2017-03-03	12	36.16
12,679	2017-03-16	2017-03-20	5	35.13
12,991	2018-02-02	2018-02-08	7	36.14
12,995	2018-02-14	2018-02-17	4	35.45
13,682	2020-02-05	2020-02-09	5	34.81
13,685	2020-02-14	2020-02-17	4	34.63
13,686	2020-02-19	2020-02-25	7	36.15
13,687	2020-02-27	2020-03-01	4	36.11
13,690	2020-03-07	2020-03-09	3	36.73
13,694	2020-03-14	2020-03-17	4	34.75
13,719	2020-04-13	2020-04-15	3	33.16
13,974	2020-12-29	2020-12-31	3	33.91
13,979	2021-01-08	2021-01-13	6	34.11
13,980	2021-01-15	2021-01-19	5	33.99
13,985	2021-01-25	2021-01-27	3	34.20
13,986	2021-01-29	2021-02-01	4	34.06
13,989	2021-02-09	2021-02-12	4	35.34
13,993	2021-02-20	2021-02-22	3	33.90
13,998	2021-02-28	2021-03-07	8	36.32
14,017	2021-03-27	2021-04-05	10	35.10
14,325	2022-02-08	2022-02-10	3	34.52
14,328	2022-02-17	2022-02-19	3	35.27
14,329	2022-02-21	2022-02-25	5	35.42
14,351	2022-03-25	2022-03-31	7	35.63
14,639	2023-01-18	2023-01-26	9	35.44
14,645	2023-02-02	2023-02-12	11	35.18
14,651	2023-02-19	2023-02-26	8	34.91
14,913	2023-11-30	2023-12-02	3	34.84
14,914	2023-12-04	2023-12-07	4	34.57
14,915	2023-12-09	2023-12-11	3	34.39
14,917	2023-12-14	2023-12-31	18	35.74

Summarize Longest Heatwave Per Year

# Extract year from start date of each heatwave
heatwave_events <- heatwave_events %>%
  mutate(Year = year(Start))

# For each year, find the longest heatwave
longest_heatwave_yearly <- heatwave_events %>%
  group_by(Year) %>%
  summarise(
    Start = Start[which.max(Duration)],
    End = End[which.max(Duration)],
    Duration = max(Duration, na.rm = TRUE),
    Max_Temp = max(Max_Temp, na.rm = TRUE),
    .groups = "drop"
  )

longest_heatwave_yearly

# A tibble: 22 × 5
    Year Start      End        Duration Max_Temp
   <dbl> <date>     <date>        <int>    <dbl>
 1  1983 1983-03-17 1983-03-20        4     34.2
 2  1987 1987-03-02 1987-03-06        5     34.9
 3  1988 1988-02-19 1988-02-22        4     34.5
 4  1991 1991-01-29 1991-02-06        9     34.1
 5  1992 1992-02-05 1992-02-23       19     38.2
 6  1998 1998-02-04 1998-02-08        5     36.0
 7  1999 1999-12-18 1999-12-23        6     34.3
 8  2000 2000-02-11 2000-03-25       44     37.6
 9  2002 2002-03-03 2002-03-06        4     34.2
10  2004 2004-03-19 2004-03-22        4     34.2
# ℹ 12 more rows

Export for Publication

# CSV
write.csv(longest_heatwave_yearly, "Nsukka_Longest_Heatwaves_Per_Year.csv", row.names = FALSE)

Create a Nice Table

flextable(longest_heatwave_yearly) %>%
  set_caption("Table: Longest Heatwave Events per Year in Nsukka") %>%
  autofit()

Year	Start	End	Duration	Max_Temp
1,983	1983-03-17	1983-03-20	4	34.24
1,987	1987-03-02	1987-03-06	5	34.87
1,988	1988-02-19	1988-02-22	4	34.49
1,991	1991-01-29	1991-02-06	9	34.06
1,992	1992-02-05	1992-02-23	19	38.16
1,998	1998-02-04	1998-02-08	5	35.97
1,999	1999-12-18	1999-12-23	6	34.33
2,000	2000-02-11	2000-03-25	44	37.57
2,002	2002-03-03	2002-03-06	4	34.19
2,004	2004-03-19	2004-03-22	4	34.17
2,005	2005-02-01	2005-02-04	4	34.39
2,006	2006-02-06	2006-02-10	5	35.15
2,007	2007-03-17	2007-03-20	4	34.33
2,010	2010-03-28	2010-03-30	3	34.23
2,015	2015-02-28	2015-03-08	9	35.66
2,016	2016-01-29	2016-03-07	39	37.87
2,017	2017-02-20	2017-03-03	12	36.16
2,018	2018-02-02	2018-02-08	7	36.14
2,020	2020-02-19	2020-02-25	7	36.73
2,021	2021-03-27	2021-04-05	10	36.32
2,022	2022-03-25	2022-03-31	7	35.63
2,023	2023-12-14	2023-12-31	18	35.74

Plot Longest Heatwave Duration Per Year

ggplot(longest_heatwave_yearly, aes(x = Year, y = Duration)) +
  geom_line(color = "red", size = 1) +
  geom_point(color = "darkred", size = 2) +
  geom_smooth(method = "lm", se = TRUE, color = "blue", linetype = "dashed") +
  labs(title = "Trend of Longest Heatwave Duration in Nsukka",
       x = "Year", y = "Duration (days)") +
  theme_minimal()

Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.

`geom_smooth()` using formula = 'y ~ x'

Plot Maximum Temperature of Longest Heatwave Per Year

ggplot(longest_heatwave_yearly, aes(x = Year, y = Max_Temp)) +
  geom_line(color = "orange", size = 1) +
  geom_point(color = "darkorange", size = 2) +
  geom_smooth(method = "lm", se = TRUE, color = "blue", linetype = "dashed") +
  labs(title = "Trend of Maximum Temperature During Longest Heatwave (Nsukka)",
       x = "Year", y = "Max Temperature (°C)") +
  theme_minimal()

`geom_smooth()` using formula = 'y ~ x'

(Optional) Statistical Test for Trend

# Duration trend
model_duration <- lm(Duration ~ Year, data = longest_heatwave_yearly)
summary(model_duration)


Call:
lm(formula = Duration ~ Year, data = longest_heatwave_yearly)

Residuals:
   Min     1Q Median     3Q    Max 
-7.974 -5.510 -3.838 -0.510 34.321 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)
(Intercept) -249.2255   392.3717  -0.635    0.533
Year           0.1295     0.1956   0.662    0.516

Residual standard error: 11.12 on 20 degrees of freedom
Multiple R-squared:  0.02142,   Adjusted R-squared:  -0.0275 
F-statistic: 0.4379 on 1 and 20 DF,  p-value: 0.5157

# Max temp trend
model_temp <- lm(Max_Temp ~ Year, data = longest_heatwave_yearly)
summary(model_temp)


Call:
lm(formula = Max_Temp ~ Year, data = longest_heatwave_yearly)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.3964 -0.9060 -0.3399  0.3014  3.1676 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)
(Intercept) -35.16752   43.62530  -0.806    0.430
Year          0.03522    0.02175   1.619    0.121

Residual standard error: 1.236 on 20 degrees of freedom
Multiple R-squared:  0.1159,    Adjusted R-squared:  0.0717 
F-statistic: 2.622 on 1 and 20 DF,  p-value: 0.121

Create the Two Plots a combined figure with side-by-side plots (Duration vs Max Temp) that will be neat and publication-ready.

# Plot 1: Duration Trend
p1 <- ggplot(longest_heatwave_yearly, aes(x = Year, y = Duration)) +
  geom_line(color = "red", size = 1) +
  geom_point(color = "darkred", size = 2) +
  geom_smooth(method = "lm", se = TRUE, color = "blue", linetype = "dashed") +
  labs(title = "Trend of Longest Heatwave Duration",
       x = "Year", y = "Duration (days)") +
  theme_minimal()

# Plot 2: Max Temp Trend
p2 <- ggplot(longest_heatwave_yearly, aes(x = Year, y = Max_Temp)) +
  geom_line(color = "orange", size = 1) +
  geom_point(color = "darkorange", size = 2) +
  geom_smooth(method = "lm", se = TRUE, color = "blue", linetype = "dashed") +
  labs(title = "Trend of Maximum Temperature During Longest Heatwave",
       x = "Year", y = "Max Temp (°C)") +
  theme_minimal()


# Step 3: Combine into One Figure
 
# Side-by-side
combined_plot <- p1 + p2

# Show plot
combined_plot

`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

# stacked (one above the other)
combined_plot <- p1 / p2
combined_plot

`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

Save for Publication

ggsave("Nsukka_Heatwave_Trends.png", combined_plot, width = 12, height = 6, dpi = 300)

`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

ggsave("Nsukka_Heatwave_Trends.pdf", combined_plot, width = 12, height = 6)

`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

Annotate the hottest year and longest duration so plots are eye-catching and publication-ready.

Highlight:

The year with the longest heatwave duration
The year with the highest max temperature

Identify Extremes

# Longest duration
longest_event <- longest_heatwave_yearly[which.max(longest_heatwave_yearly$Duration), ]

# Hottest event
hottest_event <- longest_heatwave_yearly[which.max(longest_heatwave_yearly$Max_Temp), ]

Annotated Plots

# Plot 1: Duration Trend with annotation
p1 <- ggplot(longest_heatwave_yearly, aes(x = Year, y = Duration)) +
  geom_line(color = "red", size = 1) +
  geom_point(color = "darkred", size = 2) +
  geom_smooth(method = "lm", se = TRUE, color = "blue", linetype = "dashed") +
  geom_point(data = longest_event, aes(x = Year, y = Duration),
             color = "black", size = 3, shape = 21, fill = "yellow") +
  geom_text(data = longest_event, aes(x = Year, y = Duration,
             label = paste0("Longest: ", Duration, " days (", Year, ")")),
             vjust = -1, hjust = 0.5, color = "black", size = 4) +
  labs(title = "Trend of Longest Heatwave Duration",
       x = "Year", y = "Duration (days)") +
  theme_minimal()

# Plot 2: Max Temp Trend with annotation
p2 <- ggplot(longest_heatwave_yearly, aes(x = Year, y = Max_Temp)) +
  geom_line(color = "orange", size = 1) +
  geom_point(color = "darkorange", size = 2) +
  geom_smooth(method = "lm", se = TRUE, color = "blue", linetype = "dashed") +
  geom_point(data = hottest_event, aes(x = Year, y = Max_Temp),
             color = "black", size = 3, shape = 21, fill = "yellow") +
  geom_text(data = hottest_event, aes(x = Year, y = Max_Temp,
             label = paste0("Hottest: ", Max_Temp, "°C (", Year, ")")),
             vjust = -1, hjust = 0.5, color = "black", size = 4) +
  labs(title = "Trend of Maximum Temperature During Longest Heatwave",
       x = "Year", y = "Max Temp (°C)") +
  theme_minimal()

Combine & Save

# Side by side
combined_plot <- p1 + p2

# Save high-res
ggsave("Nsukka_Heatwave_Trends_Annotated.png", combined_plot, width = 12, height = 6, dpi = 300)

`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

ggsave("Nsukka_Heatwave_Trends_Annotated.pdf", combined_plot, width = 12, height = 6)

`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

Count Heatwave Events Per Year

# Count heatwave events per year
heatwave_frequency <- heatwave_events %>%
  mutate(Year = year(Start)) %>%
  group_by(Year) %>%
  summarise(Num_Events = n(),
            Avg_Duration = mean(Duration, na.rm = TRUE),
            Max_Duration = max(Duration, na.rm = TRUE),
            .groups = "drop")

heatwave_frequency

# A tibble: 22 × 4
    Year Num_Events Avg_Duration Max_Duration
   <dbl>      <int>        <dbl>        <int>
 1  1983          1          4              4
 2  1987          4          4.5            5
 3  1988          2          4              4
 4  1991          2          6              9
 5  1992          4         11.5           19
 6  1998          4          4              5
 7  1999          3          4              6
 8  2000          8         11.1           44
 9  2002          1          4              4
10  2004          2          3.5            4
# ℹ 12 more rows

Plot Frequency Over Time

ggplot(heatwave_frequency, aes(x = Year, y = Num_Events)) +
  geom_line(color = "firebrick", size = 1) +
  geom_point(color = "darkred", size = 2) +
  geom_smooth(method = "lm", se = TRUE, color = "blue", linetype = "dashed") +
  labs(title = "Trend of Heatwave Frequency in Nsukka",
       x = "Year", y = "Number of Heatwave Events") +
  theme_minimal()

`geom_smooth()` using formula = 'y ~ x'

Statistical Trend Test

model_freq <- lm(Num_Events ~ Year, data = heatwave_frequency)
summary(model_freq)


Call:
lm(formula = Num_Events ~ Year, data = heatwave_frequency)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.1930 -1.5485 -0.1039  1.5514  4.7701 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)  
(Intercept) -189.38729   78.85706  -2.402   0.0262 *
Year           0.09631    0.03932   2.450   0.0236 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.235 on 20 degrees of freedom
Multiple R-squared:  0.2308,    Adjusted R-squared:  0.1923 
F-statistic:     6 on 1 and 20 DF,  p-value: 0.02364

Export for Publication

# Export frequency table
write.csv(heatwave_frequency, "Nsukka_Heatwave_Frequency.csv", row.names = FALSE)