library(terra)
library(sf)
library(stars)
library(ggplot2)
library(viridis)
#Step1: Read the raster
b <- terra::rast("./temp/cru_ts4.07.1901.2022.tmp.dat.nc")
#Step2: Extract time information
time_values <- time(b)
#Step3: Turning strings into time
date_objects <- as.Date(time_values)
# Step 4: Select only the dates from the year 1972 and 2022 
dates_1972 <- date_objects[format(date_objects, "%Y") == "1972"]
dates_2022 <- date_objects[format(date_objects, "%Y") == "2022"]
# Step 5: Find the index of the selected dates in the time_values
selected_indices_1972 <- which(time_values %in% dates_1972)
selected_indices_2022 <- which(time_values %in% dates_2022)
# Step 6: Extract raster values for the selected dates
selected_rasters_1972 <- b[[selected_indices_1972]]
selected_rasters_2022 <- b[[selected_indices_2022]]
# Step 7: Extract layer names
layer_names_1972 <- names(selected_rasters_1972)
layer_names_2022 <- names(selected_rasters_2022)
# Step 8: Find indices of layers that start with "tmp"
indices_tmp_1972 <- grep("^tmp", layer_names_1972)
tmp_layers_1972 <- selected_rasters_1972[[indices_tmp_1972]]
indices_tmp_2022 <- grep("^tmp", layer_names_2022)
tmp_layers_2022 <- selected_rasters_2022[[indices_tmp_2022]]
# Step 9: Changing raster names to temperature
# Ensure that the number of names matches the number of layers
if(length(indices_tmp_1972) == length(dates_1972)) {
  names(tmp_layers_1972) <- format(dates_1972, "%Y-%m-%d")
}

if(length(indices_tmp_2022) == length(dates_2022)) {
  names(tmp_layers_2022) <- format(dates_2022, "%Y-%m-%d")
}
#Step10: Turning the spatraster into a list of raster
stars_list_1972 <- lapply(tmp_layers_1972, function(raster_obj) st_as_stars(raster::raster(raster_obj)))
stars_list_2022 <- lapply(tmp_layers_2022, function(raster_obj) st_as_stars(raster::raster(raster_obj)))
#Step11: Use lapply to apply st_apply to each stars object to calculate the mean
mean_list_1972 <- lapply(stars_list_1972, function(stars_obj) st_apply(stars_obj, c("x", "y"), mean))
mean_list_2022 <- lapply(stars_list_2022, function(stars_obj) st_apply(stars_obj, c("x", "y"), mean))
#Step12: Combine the list of results into a single stars object
mean_stars_1972 <- do.call(stars::st_as_stars, mean_list_1972)
mean_stars_2022 <- do.call(stars::st_as_stars, mean_list_2022)

#Plotting
fig <-ggplot() +
  geom_stars(data = mean_stars_1972) +
  scale_fill_viridis_c( na.value = NA) +
  theme_bw() +
  ggtitle("1972") +
    labs(fill = "" )+
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    axis.title.x = element_blank(),  
    axis.title.y = element_blank())

#Plotting
fig2<-ggplot()+
  geom_stars(data=mean_stars_2022)+
  scale_fill_viridis_b(na.value=NA)+
  theme_bw()+
  ggtitle(2022)+
  labs(fill = "" )+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
    axis.title.x = element_blank(),  
    axis.title.y = element_blank())

# Reading the file
average_temp <- read.csv("./average_temp_1901_2022.csv")

# Subsetting to between 1972 and 2022 (inclusive)
average_temp_1972_2022 <- subset(average_temp, year >= 1972 & year <= 2022)

# Calculating the 1972-2022 average
average_temp_1972_2022_avg <- mean(average_temp_1972_2022$average_temperature)

# Calculate the deviations from the 1972-2022 average only for the years 1972 to 2022
average_temp_1972_2022$dev_1972_2022_avg <- average_temp_1972_2022$average_temperature - average_temp_1972_2022_avg

# Create a new column for bar color based on the condition, but only for the years 1972 to 2022
average_temp_1972_2022$bar_color <- ifelse(average_temp_1972_2022$dev_1972_2022_avg > 0, "red", "blue")

# Add the deviation and bar color columns to the original dataframe for years 1972 to 2022
average_temp <- merge(average_temp, average_temp_1972_2022[, c("year", "dev_1972_2022_avg", "bar_color")], by = "year", all.x = TRUE)

# Create the bar plot
v<-ggplot(average_temp, aes(x = year, y = dev_1972_2022_avg, fill = bar_color)) +
  geom_bar(stat = "identity", color = "black") +
  scale_fill_manual(values = c("blue", "red"), guide="none") +
  theme_bw() +
  labs(title = "Time Trend of Average Temperature",
       x = "Year",
       y = "Difference from 1972-2022 avearge (degrees C)")

library(sf)
library(ggplot2)
library(units)
library(ggpubr)
# Reading borders at level 0 and 2
nz0 <- st_read(dsn="./gadm41_NZL_shp/gadm41_NZL_0.shp", quiet = TRUE)
nz2 <- st_read(dsn="./gadm41_NZL_shp/gadm41_NZL_2.shp", quiet = TRUE)
#Simplifying the shapefile for easier plotting
nz0<-st_simplify(nz0, dTolerance=0.02)
nz2<-st_simplify(nz2, dTolerance=0.005)
#Extracting only the relevant variable
nz2b<-subset(nz2, select = c("NAME_2"))
#Calculating the area
nz2b$area<-st_area(nz2b)
nz2b$area_sqkm<-set_units(nz2b$area, 'km^2')
nz2b$area_sqkm_numeric<-as.numeric(nz2b$area_sqkm)
#setting extents
min_lon_x <- 166
max_lon_x <- 179   
min_lat_y <- - 48  
max_lat_y <- -34
#1972
#Ensuring that the CRS is the same
st_crs(mean_stars_1972)<-st_crs(nz0)
# Crop the resampled raster data to NZ boundaries
cropped_1972<- st_crop(mean_stars_1972, nz0)
# Aggregating lights by polygon
agg_1972<- raster::aggregate(cropped_1972, nz2, mean, na.rm=TRUE)

f_1972<- ggplot()+
  geom_stars(data=agg_1972)+
  coord_sf(
    xlim = c(min_lon_x , max_lon_x),  
    ylim = c(min_lat_y , max_lat_y)  
  ) +
    ggtitle("1972")+
    labs(fill = "" ) +
    scale_fill_viridis_b(na.value = NA)

#2022
# Ensuring that the CRS is the same
st_crs(mean_stars_2022)<-st_crs(nz0)
# Crop the resampled raster data to China's boundaries
cropped_2022<- st_crop(mean_stars_2022, nz0)
#Step2: Aggregating lights by polygon
agg_2022<- raster::aggregate(cropped_2022, nz2, mean, na.rm=TRUE)

f_2022<-ggplot()+
  geom_stars(data=agg_2022)+
  coord_sf(
    xlim = c(min_lon_x , max_lon_x),  
    ylim = c(min_lat_y , max_lat_y)  
  ) +
    ggtitle("2022")+
    labs(fill = "" ) +
  scale_fill_viridis_b(na.value = NA)
# Comparing the Two Plots
mean_temp<-ggarrange(f_1972, f_2022, common.legend = TRUE)
arranged_mean <- annotate_figure(mean_temp,top = text_grob("Average Temperature 1972 vs.2022", 
                                                             size = 14, 
                                                             face = "bold"))

library(ggplot2)
library(viridis)
# Calculate the difference
diff_raster <- agg_2022 - agg_1972 
# Plot the difference using ggplot and viridis theme
diff_map <- ggplot() +
  geom_stars(data = diff_raster) +
  scale_fill_viridis_c(option="turbo", name = "Temp.
  2022-1972", na.value = NA) +
  theme_bw()+
  coord_sf(
    xlim = c(min_lon_x , max_lon_x),  
    ylim = c(min_lat_y , max_lat_y)  
  ) +
theme(legend.position = c(0, 1),  
        legend.justification = c(0, 1), 
        legend.background = element_rect(color = "black", size = 0.2))  

library(sf)
sf_use_s2(FALSE)
library(ggplot2)
library(ggrepel)
# Loading New Zealand Shapefile
nz_cntry1 <- st_read(dsn="./gadm41_NZL_shp/gadm41_NZL_1.shp", quiet = TRUE)
# Simplify lines
nz1<-st_simplify(nz_cntry1,  dTolerance = 0.05)
#Extracting only the relevant variable
nz1b<-subset(nz1, select = c("NAME_1"))
#Calculating the area
nz1b$area<-st_area(nz1b)
nz1b$area_sqkm<-set_units(nz1b$area, 'km^2')
nz1b$area_sqkm_numeric<-as.numeric(nz1b$area_sqkm)

library(ggrepel)
# Extract centroids
centroids <- st_centroid(nz1b)
# Convert to a simple data frame for ggplot
centroids_df <- as.data.frame(st_coordinates(centroids))
centroids_df$NAME_1 <- nz1b$NAME_1
centroids_df$area_sqkm_numeric <- nz1b$area_sqkm_numeric

k<-ggplot(nz1b, aes(fill = area_sqkm_numeric)) + 
  geom_sf() + 
  geom_label_repel(data = centroids_df, aes(x = X, y = Y, label = NAME_1), 
                   colour = "black", size = 2, fill = "white",
                   box.padding = unit(0.6, "lines")) +  # Adjust box.padding as necessary
  coord_sf(xlim = c(min_lon_x, max_lon_x), ylim = c(min_lat_y, max_lat_y)) +
  scale_fill_viridis_c(option = "viridis")+
  theme_void()

library(dplyr)
#Step1: Reading droughts
droughts <- read.csv('./drought.data/drought-spei-1962-2022.csv')
#Step2: Turning the CSV file into an SF object
droughts_sf<-st_as_sf(droughts, coords = c("lon", "lat"))
#Step4: Associating the crs information for the NZ shapefile to points
droughts_sf <- st_set_crs(droughts_sf, st_crs(nz1))  
#1972
droughts_1972 <- subset(droughts_sf, year == "1972")
#2022
droughts_2022 <- subset(droughts_sf, year == "2022")

library(sf)
sf_use_s2(FALSE)
library(ggplot2)

plotting_nz<-ggplot() +
  geom_sf(data = nz1b, fill = "white", color = "black", size = NA) +  # Base map
  geom_sf(data = droughts_1972, aes(size = intensity, color = intensity), alpha = 0.7) +  # Drought data
  scale_color_gradient2(
    low = "blue", mid = "yellow", high = "red",  
    breaks = c(0.25, 0.50, 0.75),
    name = "Intensity"
  ) +
  scale_size_continuous(
    breaks = c(0.25, 0.50, 0.75),  # Size breaks matching color breaks
    range = c(0.2, 3),
    name = "Intensity"
  ) +
  coord_sf(
    xlim = c(min_lon_x, max_lon_x),  
    ylim = c(min_lat_y, max_lat_y)  
  ) +
  theme_minimal() +
  labs(x = "Longitude", y = "Latitude", title = "1972: Drought Intensity in New Zealand") +
  guides(
    color = guide_legend(),
    size = guide_legend()
  )

library(sf)
sf_use_s2(FALSE)
library(ggplot2)


plotting_22 <- ggplot() +
  geom_sf(data = nz1b, fill = "white", color = "black", size = NA) +  # Base map
  geom_sf(data = droughts_2022, aes(size = intensity, color = intensity), alpha = 0.7) +  # Drought data
  scale_color_gradient2(
    low = "blue", mid = "yellow", high = "red",  # Adjusted midpoint and colors
    breaks = c(0.25, 0.50, 0.75),
    name = "Intensity"
  ) +
  scale_size_continuous(
    breaks = c(0.25, 0.50, 0.75), 
    range = c(0.2, 3),
    name = "Intensity"
  ) +
  coord_sf(
    xlim = c(min_lon_x, max_lon_x),  
    ylim = c(min_lat_y, max_lat_y)  
  ) +
  theme_minimal() +
  labs(x = "Longitude", y = "Latitude", title = "2022: Drought Intensity in New Zealand") +
  guides(
    color = guide_legend(),
    size = guide_legend()
  )

# Comparing the Two Results
arranged_plots <- ggarrange(fig3, fig4, common.legend = TRUE)
arranged_plots_title <- annotate_figure(arranged_plots, 
                                             top = text_grob("Drought Intensity per Region 1972 vs.2022", 
                                                             size = 14, 
                                                           face = "bold"))

library(sf)

# Simplify geometries
droughts_sf <- st_simplify(droughts_sf, preserveTopology = TRUE, dTolerance = 0.5)

# Filter data 
droughts_sf <- droughts_sf %>%
  filter(year >= 1972)  

# Load required libraries
library(ggplot2)
library(gganimate)
library(sf)
library(dplyr)
library(transformr)  

p <- ggplot() +
  geom_sf(data = nz1b, fill = "white", color = "black") +
  geom_sf(data = droughts_sf, aes(size=intensity, color = intensity), alpha = 0.7) +
  scale_color_gradient2(
    low = "blue", mid = "yellow", high = "red", 
    breaks = c(0.25, 0.50, 0.75),
    name = "Intensity"
  ) +
  scale_size_continuous(
    breaks = c(0.25, 0.50, 0.75),  
    range = c(0.2, 3),
    name = "Intensity"
  ) +
  coord_sf(
    xlim = c(min_lon_x, max_lon_x),  
    ylim = c(min_lat_y, max_lat_y)  
  ) +
  theme_minimal() +
  labs(x = "Longitude", y = "Latitude", title = " Drought Intensity in New Zealand") +
  guides(
    color = guide_legend(),
    size = guide_legend()
  )+
  labs(title = 'Year: {current_frame}') +
  transition_manual(year)

#1972
x1<-st_join(droughts_1972,nz1b)
xdf2<- st_drop_geometry(x1)
sdf3<-xdf2%>%
  group_by(NAME_1)%>%
  summarize(mean_int = mean(intensity, na.rm=TRUE))

fin_df72<-left_join(nz1b, sdf3, by = c('NAME_1'='NAME_1'))

#Plotting
fig3 <- ggplot() +
  geom_sf(data = fin_df72, aes(fill = mean_int), linewidth = 0.07) +
  scale_fill_viridis_b(na.value = NA, name = "mean intensity") +  # Update the name for the legend
  labs(
    title = "1972",
    fill = "Drought Intensity",  # Optionally adjust this if you want a different description
    x = "",
    y = ""
  ) +
  coord_sf(
    xlim = c(min_lon_x, max_lon_x),  
    ylim = c(min_lat_y, max_lat_y)  
  ) +
  theme_bw()

#2022
x<-st_join(droughts_2022,nz1b)

xdf<- st_drop_geometry(x)
sdf2<-xdf%>%
  group_by(NAME_1)%>%
  summarize(mean_int_2 = mean(intensity, na.rm=TRUE))

fin_df<-left_join(nz1b, sdf2, by = c('NAME_1'='NAME_1'))

#Plotting
fig4<-ggplot() +
  geom_sf(data = fin_df, aes(fill = mean_int_2), linewidth = 0.07) +
  scale_fill_viridis_b(name = "mean intensity", na.value = NA) +  
  labs(title = "2022", fill = "Drought Intensity") +  
  coord_sf(
    xlim = c(min_lon_x , max_lon_x),  
    ylim = c(min_lat_y , max_lat_y)  
  ) +
  labs(x = "", y = "")+
  theme_bw()

# Comparing the Two Results
arranged_plots <- ggarrange(fig3, fig4, common.legend = TRUE)
arranged_plots_title <- annotate_figure(arranged_plots, 
                                             top = text_grob("Drought Intensity per Region 1972 vs.2022", 
                                                             size = 14, 
                                                           face = "bold"))

library(sf)
library(dplyr)
library(ggplot2)
library(gganimate)
 

# Joining the drought data with the spatial data for New Zealand regions
combined_data <- st_join(droughts_sf, nz1b)

# Dropping geometry for the calculation, then grouping and summarizing by NAME_1 and year
mean_intensity_by_year <- combined_data %>%
  group_by(NAME_1, year) %>%
  summarize(mean_int = mean(intensity, na.rm = TRUE), .groups = 'drop')

# Joining summarized data back to the spatial data
animated_data <- st_join(nz1b, mean_intensity_by_year, by = 'NAME_1')
z_anim <- ggplot() +
  geom_sf(data = animated_data, aes(fill = mean_int), size = 0.1) +
  scale_fill_viridis_b(option="viridis", name = "Drought Intensity",  na.value = NA) +
  labs(title = "Year: {frame_time}", x = "", y = "") +
  theme_minimal() +
  coord_sf(
    xlim = c(min_lon_x , max_lon_x),  
    ylim = c(min_lat_y , max_lat_y)  
  ) +
  labs(title = 'Year: {current_frame}') +
  transition_manual(year)