---
title: "Why does Australia need to transition to renewable energy?"
author: "Benjamin Johnson S3930970"
date: "`r Sys.Date()`"
output: 
  flexdashboard::flex_dashboard:
    orientation: columns
    vertical_layout: fill
    # storyboard: true
    theme: cerulean
    social: ["twitter", "facebook", "linkedin"]
    source_code: embed
always_allow_html: yes
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)

#Question: "Why does Australia need to transition to renewable energy?"

#install.packages("ncdf4")
#install.packages("raster")
#install.packages("rgdal")
#install.packages("tidyverse")
#install.packages("chron")
#install.packages("flexdashboard")
#install.packages("leaflet")
#install.packages("geojsonio")
#install.packages("htmltools")
#install.packages("fontawesome")
#install.packages("htmlwidgets")

# if(!require('plotly')) {
#   install.packages('plotly')
#   library('plotly')
# }

library(ggplot2)
library(dplyr)
library(knitr)
library(readr) # Useful for importing data
library(magrittr) # Enables the use of pipe operators for readability
library(tidyr) # Reshapes and splits data/cells
library(stringr) # Handy string operations
library(ggthemes) #Additional ggplot themes
library(plotly) #Enable plotly for interaction
# library(ncdf4) # package for netcdf manipulation
library(raster) # package for raster manipulation
# library(rgdal) # package for geospatial analysis
library(tidyverse) # tidyverse things
library(chron) #for time conversion
library(flexdashboard) #For use of Flex dashboard
library(leaflet) #Easier for mapping
library(geojsonio) #For GeoJSON use
library(htmltools) #For formatting
library(htmlwidgets) #For formatting
library(countrycode) #To look up the ISO code for each country
library(openxlsx) #To import xlsx
#library(fontawesome) #To use additional icons

```

1. Introduction
================================================================

Column {data-width="1%"}
-----------------------------------------------------------------------

Column {data-width="70%"}
-----------------------------------------------------------------------
<br>
This year, the latest report from the CSIRO (Australia's National Science Agency) have noted a few key points from Dr. Karl Braganza [State of the Climate Report 2022]:

* "Australia's climate has warmed on average by **1.47** degrees since 1910..." (Degrees celsius which will be the consistent measure used throughout)
* "...contrasting rainfall trends across the north and the south of the country." with one of the most significant flood periods during La Niña in 2021-22.
* “We’re expecting to see longer fire seasons in the future for the south and east, and an increase in the number of dangerous fire weather days..."

All of which spells an ominous future but *how did we get here?* 
And *why is there so much focus on "clean energy"?*

Through the following storyboard we look to journey through:

1. Global Temperatures (and the trend through the years) to... 
2. What is causing changes in temperatures and from there... 
3. Looking for the biggest contributor to this change in Australia.
4. And lastly how we might be able to help reduce any impact.
<br>
<br>

**Click on the 'Analysis via Storyboard' tab at the top to continue.**
<br>
<br>

***

**Reference:** CSIRO (Australia’s National Science Agency), *‘What are the impacts of extreme weather and climate events?’*. Accessed 24 November 2022, CSIRO website: https://www.csiro.au/en/news/News-releases/2022/State-of-the-Climate-Report-2022

Column {data-width="3%"}
-----------------------------------------------------------------------

Column {data-width="25%"}
-----------------------------------------------------------------------

```{r fig.cap = "Source: Lake Eildon - by Benjamin Johnson", echo=FALSE} 

knitr::include_graphics("_DSC3955.jpg")

```

Column {data-width="1%"}
-----------------------------------------------------------------------

2. Analysis via Storyboard {.storyboard}
================================================================

### How are temperatures now compared to 70 years ago? <br>(1/8)

```{r fig.width = 21, fig.height = 7, fig.align = "center", echo=FALSE}

# {data-commentary-width=350}

GlobalTempAnomaly <- read_csv("latest.txt", col_names = FALSE)

#<start year><start month><end year><end month>   <mean temperature anomaly (°C)>
#Need to split out the single variable as per the above

GlobalTempAnomaly %<>% 
  separate("X1", into = c("date", "mean_temp_anomaly"), sep = "   ")

GlobalTempAnomaly %<>% 
  mutate(start_year = as.numeric(substr(GlobalTempAnomaly$date, 1, 4)),
         start_month = as.numeric(substr(GlobalTempAnomaly$date, 5, 6)),
         end_year = as.numeric(substr(GlobalTempAnomaly$date, 7, 10)),
         end_month = as.numeric(substr(GlobalTempAnomaly$date, 11, 12)))

GlobalTempAnomaly$mean_temp_anomaly <- as.numeric(GlobalTempAnomaly$mean_temp_anomaly)

#Replicating the same figure "derived"
GlobalTempAnomaly %<>%
  mutate(absolute_temp = mean_temp_anomaly+13.97)

#Filtering out a range prior to this storyboard
#Could display a longer period but might disconnect viewers with constant changes in time series
#Also introducing a long period with a restricted ratio might over-exaggerate any increase over time
GlobalTempAnomaly <- GlobalTempAnomaly[GlobalTempAnomaly$end_year >= 1950,]

#Validating that the derived temp average is equivalent to the baseline figure
#mean(GlobalTempAnomaly[GlobalTempAnomaly$end_year < 1991 & GlobalTempAnomaly$end_year > 1960,]$derived_mean_temp)

#summary(GlobalTempAnomaly)

# p1 <- ggplot(data = GlobalTempAnomaly[GlobalTempAnomaly$end_year < 2007,], aes(x=end_year, y=mean_temp_anomaly, colour = mean_temp_anomaly))
#   p1 + geom_line(stat = "identity")+
#   labs(x='Year', y='Mean Temperature Anomaly (°C)', title='Annual Mean Temperature Anomalies from 1850 to 2006 (Global)', subtitle = 'Source: Bureau of Meteorology - Global climate variability & change\nTime series data using 1961-1990 baseline average of 13.97 °C') +
#   theme_stata() +
#   theme(plot.title = element_text(size = 18, face="bold"),
#         plot.subtitle = element_text(size = 16, face="bold.italic"),
#         axis.text=element_text(size=15),
#         axis.title=element_text(size=16,face="bold"),
#         axis.text.y = element_text(angle = 0),
#         axis.title.y = element_text(vjust = 5),
#         panel.grid = element_line(size = 1)) +
#     scale_x_continuous(breaks = round(seq(min(GlobalTempAnomaly$end_year), max(GlobalTempAnomaly$end_year), by = 25),1)) +
#     theme(legend.position="right") +
#     scale_colour_gradient(name="Mean Temperature Anomaly (°C)\n(Baseline 1961-1990)",low="blue",high="red") +
#     geom_smooth(colour = "black", method = "lm") +
#     geom_hline(aes(yintercept = 0), size = 2) +
#     ylim(-0.6,0.6)
    
#Previous bar version but introduces some inconsistency in terms of view  
# p1 <- plot_ly(data = GlobalTempAnomaly, x = ~end_year, y = ~mean_temp_anomaly, type = "bar", color = ~mean_temp_anomaly,
#               colors = c("#67a9cf","#ef8a62"),
#               text = paste("<b>Year</b> = ", GlobalTempAnomaly$end_year,
#                            "<br><b>Mean Temp Anomaly difference (°C)</b> = ", GlobalTempAnomaly$mean_temp_anomaly,
#                            "<br><b>Absolute Value (°C)</b> = ", GlobalTempAnomaly$absolute_temp)) %>%
#   layout(title ="Annual Mean Temperature Anomalies from 1850 to 2019 (Global)",
#          yaxis = list(zeroline = TRUE, title = "Mean Temperature Anomaly (°C)"),
#          xaxis = list(zeroline = FALSE, title = "Year"),
#          showlegend = FALSE)
# 
# p1

#Range introduced as the average context we're working with is 13.97, so +1 and -1 in terms of y axis limit should even out the view
#Removed exactly 13 as a value to avoid overlap with 1950 x axis label

# labels <- "\u2103"
# print(labels)

m <- lm(absolute_temp~end_year,data = GlobalTempAnomaly)
    
GlobalTemp_p1 <- plot_ly(data = GlobalTempAnomaly, 
                         x = ~end_year, 
                         y = ~absolute_temp, 
                         type = "scatter", 
                         mode = "line", 
                         colors = c('#a6cee3','#b2df8a','#1f78b4'), #Not absolutely required but...
                         text = paste("<b>Year</b>: ", GlobalTempAnomaly$end_year,
                                      "<br><b>Mean Temp Anomaly difference from 1951-1980 (\u2103)</b>: ", GlobalTempAnomaly$mean_temp_anomaly,
                                      "<br><b>Absolute Value (\u2103)</b>: ", GlobalTempAnomaly$absolute_temp)) %>%
  layout(title =list(text = "<b>Average global temperatures increasing over time (1950 to 2019)</b>", font = list(size = 19), x = 0.01, yanchor =  'top'),
         yaxis = list(zeroline = TRUE, showgrid = TRUE, nticks = 10, title = "<b>Global Mean Temperature (\u2103)</b>", range = c(13.01,14.99)),
         xaxis = list(zeroline = FALSE, showgrid = FALSE, title = "<b>Year</b>"),
         showlegend = FALSE) %>%
  hide_colorbar() %>%
  add_lines(x=~end_year, y=predict(m), name="Linear", color='black')

GlobalTemp_p1

#library(colorblindr)
#cvd_grid(p2)
    
#Additional note re: absolute values

# In reply please quote: CAS-40694-W4V7Q7
# 
# Date: 11/15/2022
# 
# Dear Benjamin,
# 
# Thank you for your enquiry.
# 
# That is correct. The Australian climate variability & change timeseries graphs have always used temperature anomalies, not absolute values.
#  
# The base period is specified, and the average over that base period is provided with the graphs, so users can calculate corresponding absolute values, or rebase the dataset to a different base period, if they wish or that suits their use better.
# 
# Regards,
# 
# Climate (EP)
# Environmental Prediction Services - Climate
# Bureau of Meteorology

```
> **Data Source:** Commonwealth of Australia, Bureau of Meteorology. (2022). *Climate change – trends and extremes*. Accessed 8 November, 2022, Bureau of Meteorology website: http://www.bom.gov.au/climate/change/index.shtml#tabs=Tracker&tracker=global-timeseries&tQ=graph%3Dglobal_t%26region%3Dglobal%26season%3D0112 <br>
> **Data Information Reference:** Commonwealth of Australia, Bureau of Meteorology. (2022). *About the global temperature timeseries graphs*. Accessed 8 November, 2022, Bureau of Meteorology website: http://www.bom.gov.au/climate/change/about/global_temp_timeseries.shtml <br>
> **Colour scheme reference:**: Brewer c, Harrower M, and The Pennsylvania State University (2021). *ColorBrewer: Color Advice for Maps.*. Accessed 3 December, 2022, Published by Color Brewer 2 website: https://colorbrewer2.org/#type=qualitative&scheme=Paired&n=3

***

* Taking the 1950s onwards, we can see a steady increase in temperature despite the expected fluctuations.
* Many climate change graphs use a base period to measure how much temperature has changed. In this case we will use 1951-1980 where necessary.
* Hovering over each data point, we can also see how much the temperature has changed since that period.

**Note:** This data series extends back to 1900 but has been limited to 1950 to keep a more consistent time frame through out the storyboard. This also means the linear regression (line) is also affected having a slightly steeper slope than compared to a time period with lower temperature values (before 1950).

### That temperature increase may not necessarily impact my country though?  <br>(2/8)

```{r fig.width = 21, fig.height = 21, fig.align = "center", echo=FALSE}

#Base period 1951 - 1980

TempAnomCountry <- read_csv("Tdata.csv", col_names = TRUE)

TempAnomCountry %<>%
  pivot_longer(!c(1:3),names_to="Year", values_to = "Temp")

colnames(TempAnomCountry)[2] <- 'iso_a3'

TempAnomCountry$Year <- as.numeric(TempAnomCountry$Year)

geoData <- geojsonio::geojson_read("custom.geo.json", what = "sp")
#geoData <- readLines("custom.geo.json")
# topoData <- readLines("countries-110m.json")

#names(geoData)
#Need to use "iso_a3"

# leaflet() %>% setView(lng = 0.0, lat = 0.0, zoom = 1.5) %>%
#   addTiles() %>%
#   addTopoJSON(topoData, weight = 1, color = "#444444", fill = FALSE)

# globe_map <- leaflet() %>% setView(lng = 0.0, lat = 0.0, zoom = 1.5) %>%
#   addTiles() %>%
#   addGeoJSON(geoData, weight = 1, color = "#444444", fill = TRUE, opacity = 1, fillOpacity = 0.8)
# 
# globe_map

# globe_temps <-sp::merge(geoData, TempAnomCountry[TempAnomCountry$Year == 2019,], by="iso_a3", duplicateGeoms = TRUE)

temps_2019 <-sp::merge(geoData, TempAnomCountry[TempAnomCountry$Year == 2019,], by="iso_a3", duplicateGeoms = TRUE)
#Create the hover labels for each year...
labels_2019 <- sprintf(
  "<strong>%s %s</strong><br/>%.2g\u2103 (Mean Temperature Difference compared to the 1951-1980 mean)",
  temps_2019$Country, 
  temps_2019$Year,
  temps_2019$Temp) %>% 
  lapply(htmltools::HTML)
temps_2010 <-sp::merge(geoData, TempAnomCountry[TempAnomCountry$Year == 2010,], by="iso_a3", duplicateGeoms = TRUE)
labels_2010 <- sprintf(
  "<strong>%s %s</strong><br/>%.2g\u2103 (Mean Temperature Difference compared to the 1951-1980 mean)",
  temps_2010$Country, 
  temps_2010$Year,
  temps_2010$Temp) %>% 
  lapply(htmltools::HTML)
temps_2000 <-sp::merge(geoData, TempAnomCountry[TempAnomCountry$Year == 2000,], by="iso_a3", duplicateGeoms = TRUE)
labels_2000 <- sprintf(
  "<strong>%s %s</strong><br/>%.2g\u2103 (Mean Temperature Difference compared to the 1951-1980 mean)",
  temps_2000$Country, 
  temps_2000$Year,
  temps_2000$Temp) %>% 
  lapply(htmltools::HTML)
temps_1990 <-sp::merge(geoData, TempAnomCountry[TempAnomCountry$Year == 1990,], by="iso_a3", duplicateGeoms = TRUE)
labels_1990 <- sprintf(
  "<strong>%s %s</strong><br/>%.2g\u2103 (Mean Temperature Difference compared to the 1951-1980 mean)",
  temps_1990$Country, 
  temps_1990$Year,
  temps_1990$Temp) %>% 
  lapply(htmltools::HTML)
temps_1980 <-sp::merge(geoData, TempAnomCountry[TempAnomCountry$Year == 1980,], by="iso_a3", duplicateGeoms = TRUE)
labels_1980 <- sprintf(
  "<strong>%s %s</strong><br/>%.2g\u2103 (Mean Temperature Difference compared to the 1951-1980 mean)",
  temps_1980$Country, 
  temps_1980$Year,
  temps_1980$Temp) %>% 
  lapply(htmltools::HTML)
temps_1970 <-sp::merge(geoData, TempAnomCountry[TempAnomCountry$Year == 1970,], by="iso_a3", duplicateGeoms = TRUE)
labels_1970 <- sprintf(
  "<strong>%s %s</strong><br/>%.2g\u2103 (Mean Temperature Difference compared to the 1951-1980 mean)",
  temps_1970$Country, 
  temps_1970$Year,
  temps_1970$Temp) %>% 
  lapply(htmltools::HTML)
temps_1960 <-sp::merge(geoData, TempAnomCountry[TempAnomCountry$Year == 1960,], by="iso_a3", duplicateGeoms = TRUE)
labels_1960 <- sprintf(
  "<strong>%s %s</strong><br/>%.2g\u2103 (Mean Temperature Difference compared to the 1951-1980 mean)",
  temps_1960$Country, 
  temps_1960$Year,
  temps_1960$Temp) %>% 
  lapply(htmltools::HTML)
temps_1950 <-sp::merge(geoData, TempAnomCountry[TempAnomCountry$Year == 1950,], by="iso_a3", duplicateGeoms = TRUE)
labels_1950 <- sprintf(
  "<strong>%s %s</strong><br/>%.2g\u2103 (Mean Temperature Difference compared to the 1951-1980 mean)",
  temps_1950$Country, 
  temps_1950$Year,
  temps_1950$Temp) %>% 
  lapply(htmltools::HTML)

#Continuous
# pal <- colorNumeric("RdYlBu", domain = c(-3,-2.5,-2,-1.5,-1,-0.5,0,0.5,1,1.5,2,2.5,3), reverse = TRUE)
pal <- colorNumeric(c('#ca0020','#f4a582','#f7f7f7','#92c5de','#0571b0'), domain = c(-3,-2.5,-2,-1.5,-1,-0.5,0,0.5,1,1.5,2,2.5,3), reverse = TRUE) #Colour blind friendly palette from https://colorbrewer2.org/#type=diverging&scheme=RdBu&n=5
#Bins - reversed as Blue was used for "hot" temps
# pal <- colorBin("RdYlBu", domain = TempAnomCountry$Temp, bins = c(-3,-2,-1,0,1,2,3), reverse = TRUE)


# yearInput <- 2019
# 
# globe_p1 <- leaflet(globe_temps@data[globe_temps@data$Year == yearInput,]) %>%
#   setView(lng = 0.0, lat = 0.0, zoom = 1.5)
#Something with year selection is broken... possible reference https://stackoverflow.com/questions/51691788/nas-not-permitted-in-row-index-in-shiny-app-leaflet-map/51715088

# globe_p1 <- leaflet(globe_temps, options=leafletOptions(worldCopyJump = TRUE, minZoom = 1.5, zoomSnap = 0.1)) %>%
#   setView(0, 45, zoom=1.5) %>% 
#   addTiles()

#Create the hover labels
# labels <- sprintf(
#   "<strong>%s %s</strong><br/>%.2g°C (Mean Temperature Difference compared to 1951-1980)",
#   globe_temps$Country, 
#   globe_temps$Year,
#   globe_temps$Temp) %>% 
#   lapply(htmltools::HTML)
# 
# globe_p1 %>% 
#   addPolygons(
#   fillColor = ~pal(Temp),
#   weight = 0.7,
#   opacity = 1,
#   color = "black",
#   dashArray = "3",
#   fillOpacity = 0.5,
#   smoothFactor = 0.2,
#   stroke = FALSE,
#   highlightOptions = highlightOptions(
#     weight = 3,
#     color = "#666",
#     dashArray = "",
#     fillOpacity = 0.5,
#     bringToFront = TRUE),
#   label = labels,
#   labelOptions = labelOptions(
#     style = list("font-weight" = "normal", padding = "3px 8px"),
#     textsize = "15px",
#     direction = "auto"),
#   group = "2019") %>%
#   addLegend(pal = pal, 
#             values = ~Temp, 
#             opacity = 0.7, 
#             title = "Mean Temperature Anomaly",
#             labFormat = labelFormat(suffix = "°C"),
#             position = "bottomleft")

#Layer option experimentation

title <- tags$div(HTML("<b>Exploring global temperature increases compared to the 1951-1980 period</b><br>Adjust the year on the right side to compare temperature changes through the years."))

#Unfortunately cannot position the title in the centre? And aligning to top left conflicts with the zoom control the position of which I couldn't adjust
base_map <- leaflet(geoData, options=leafletOptions(worldCopyJump = TRUE, minZoom = 1.5, zoomSnap = 0.25)) %>%
  addControl(title, position = "topright") %>%
  setView(0, 45, zoom=2) %>% 
  addTiles()

base_map %>% 
  addPolygons(data = temps_2019,
  fillColor = ~pal(temps_2019$Temp),
  weight = 0.7,
  opacity = 1,
  color = "black",
  dashArray = "3",
  fillOpacity = 0.5,
  smoothFactor = 0.2,
  stroke = FALSE,
  highlightOptions = highlightOptions(
    weight = 3,
    color = "#666",
    dashArray = "",
    fillOpacity = 0.5,
    bringToFront = TRUE),
  label = labels_2019,
  labelOptions = labelOptions(
    style = list("font-weight" = "normal", padding = "3px 8px"),
    textsize = "12px",
    direction = "auto"),
  group = "2019") %>%
  # Add year 2010 layer
  addPolygons(data = temps_2010,
  fillColor = ~pal(temps_2010$Temp),
  weight = 0.7,
  opacity = 1,
  color = "black",
  dashArray = "3",
  fillOpacity = 0.5,
  smoothFactor = 0.2,
  stroke = FALSE,
  highlightOptions = highlightOptions(
    weight = 3,
    color = "#666",
    dashArray = "",
    fillOpacity = 0.5,
    bringToFront = TRUE),
  label = labels_2010,
  labelOptions = labelOptions(
    style = list("font-weight" = "normal", padding = "3px 8px"),
    textsize = "12px",
    direction = "auto"),
  group = "2010") %>%
  # Add year 2000 layer
  addPolygons(data = temps_2000,
  fillColor = ~pal(temps_2000$Temp),
  weight = 0.7,
  opacity = 1,
  color = "black",
  dashArray = "3",
  fillOpacity = 0.5,
  smoothFactor = 0.2,
  stroke = FALSE,
  highlightOptions = highlightOptions(
    weight = 3,
    color = "#666",
    dashArray = "",
    fillOpacity = 0.5,
    bringToFront = TRUE),
  label = labels_2000,
  labelOptions = labelOptions(
    style = list("font-weight" = "normal", padding = "3px 8px"),
    textsize = "12px",
    direction = "auto"),
  group = "2000") %>%
  # Add year 1990 layer
  addPolygons(data = temps_1990,
  fillColor = ~pal(temps_1990$Temp),
  weight = 0.7,
  opacity = 1,
  color = "black",
  dashArray = "3",
  fillOpacity = 0.5,
  smoothFactor = 0.2,
  stroke = FALSE,
  highlightOptions = highlightOptions(
    weight = 3,
    color = "#666",
    dashArray = "",
    fillOpacity = 0.5,
    bringToFront = TRUE),
  label = labels_1990,
  labelOptions = labelOptions(
    style = list("font-weight" = "normal", padding = "3px 8px"),
    textsize = "12px",
    direction = "auto"),
  group = "1990") %>%
  # Add year 1980 layer
  addPolygons(data = temps_1980,
  fillColor = ~pal(temps_1980$Temp),
  weight = 0.7,
  opacity = 1,
  color = "black",
  dashArray = "3",
  fillOpacity = 0.5,
  smoothFactor = 0.2,
  stroke = FALSE,
  highlightOptions = highlightOptions(
    weight = 3,
    color = "#666",
    dashArray = "",
    fillOpacity = 0.5,
    bringToFront = TRUE),
  label = labels_1980,
  labelOptions = labelOptions(
    style = list("font-weight" = "normal", padding = "3px 8px"),
    textsize = "12px",
    direction = "auto"),
  group = "1980") %>%
  # Add year 1970 layer
  addPolygons(data = temps_1970,
  fillColor = ~pal(temps_1970$Temp),
  weight = 0.7,
  opacity = 1,
  color = "black",
  dashArray = "3",
  fillOpacity = 0.5,
  smoothFactor = 0.2,
  stroke = FALSE,
  highlightOptions = highlightOptions(
    weight = 3,
    color = "#666",
    dashArray = "",
    fillOpacity = 0.5,
    bringToFront = TRUE),
  label = labels_1970,
  labelOptions = labelOptions(
    style = list("font-weight" = "normal", padding = "3px 8px"),
    textsize = "12px",
    direction = "auto"),
  group = "1970") %>%
  # Add year 1960 layer
  addPolygons(data = temps_1960,
  fillColor = ~pal(temps_1960$Temp),
  weight = 0.7,
  opacity = 1,
  color = "black",
  dashArray = "3",
  fillOpacity = 0.5,
  smoothFactor = 0.2,
  stroke = FALSE,
  highlightOptions = highlightOptions(
    weight = 3,
    color = "#666",
    dashArray = "",
    fillOpacity = 0.5,
    bringToFront = TRUE),
  label = labels_1960,
  labelOptions = labelOptions(
    style = list("font-weight" = "normal", padding = "3px 8px"),
    textsize = "12px",
    direction = "auto"),
  group = "1960") %>%
  # Add year 1950 layer
  addPolygons(data = temps_1950,
  fillColor = ~pal(temps_1950$Temp),
  weight = 0.7,
  opacity = 1,
  color = "black",
  dashArray = "3",
  fillOpacity = 0.5,
  smoothFactor = 0.2,
  stroke = FALSE,
  highlightOptions = highlightOptions(
    weight = 3,
    color = "#666",
    dashArray = "",
    fillOpacity = 0.5,
    bringToFront = TRUE),
  label = labels_1950,
  labelOptions = labelOptions(
    style = list("font-weight" = "normal", padding = "3px 8px"),
    textsize = "12px",
    direction = "auto"),
  group = "1950") %>%
  addLegend(pal = pal, 
            values = ~na.omit(TempAnomCountry$Temp), 
            opacity = 0.7, 
            title = "Mean Temperature Anomaly",
            labFormat = labelFormat(suffix = '\u2103'),
            position = "bottomleft") %>%
  # Layers control
  addLayersControl(
    baseGroups = c("2019", "2010","2000","1990","1980","1970","1960","1950"),
    options = layersControlOptions(collapsed = FALSE, position = "topright"))

```
>	**Transformed Data Source:** Lipponen, A (14 January 2022), *‘#Temperature anomalies 1880-2021 by country…’* [Twitter post], Antti Lipponen, accessed 29th October 2022. https://gist.github.com/anttilipp/fc2c4eb2224ff0151d2a5ac1124b585e <br>
>	**Original Data Source:**	Goddard Institute for Space Studies - NASA (National Aeronautics and Space Administration) (2021), GISTEMP Land-Ocean Temperature Index (LOTI), ERSSTv5, 1200km smoothing [data set], NASA website, accessed date by visualization author unknown. https://data.giss.nasa.gov/gistemp/ <br>
> **GeoJSON Map Source**: ashkyd (2016), GeoJSON - Download vector maps, ashkyd, accessed 21 November 2022. https://geojson-maps.ash.ms/ <br>
> **Colour scheme reference:**: Brewer c, Harrower M, and The Pennsylvania State University (2021). *ColorBrewer: Color Advice for Maps.*. Accessed 3 December, 2022, Published by Color Brewer 2 website: https://colorbrewer2.org/#type=diverging&scheme=RdBu&n=5

***

The data behind the world map view uses the average temperature from 1951 to 1980 as a grounding point, then finds out since then - how much have temperatures increased?

* Use the top right control to change the years
* Hovering over your home country, how much does the temperature change over time?
* And what is the difference now?

**Note:** The years listed in the control are only points in time - additional fluctuations could occur in the years in-between (eg. 2001 to 2009).


### What is causing these temperatures to rise?  <br>(3/8)

Referring back to the CSIRO (Australia's National Science Agency) [Why is the world warming? 2022]:

* "Human activities are increasing the concentrations of greenhouse gases in the atmosphere, leading to an 'enhanced' greenhouse effect and causing surface temperatures to rise."
* "Without any greenhouse gases, the Earth’s average surface temperature would be much lower, about -18 °C rather than today’s average of about 14 °C."
* "Concentrations of greenhouse gases continue to rise, with the global average carbon dioxide concentration now above 410 parts per million (ppm)."

<br>

**Reference:** CSIRO (Australia’s National Science Agency), *‘Why is the world warming?’*. Accessed 24 November 2022, CSIRO website: https://www.csiro.au/en/research/environmental-impacts/climate-change/Climate-change-QA/Why <br>
**Additional Reference:** Shaftel, H (2022). *Global Warming vs. Climate Change | Facts – Climate Change: Vital Signs of the Planet*. NASA website, accessed 30 October 2022. http://climate.nasa.gov/global-warming-vs-climate-change/

### Can't all these emissions just be absorbed?  <br>(4/8)

```{r fig.width = 21, fig.height = 7, fig.align = "center", echo=FALSE}

PRIMAPHist <- read_csv("Guetschow-et-al-2022-PRIMAP-hist_v2.4_11-Oct-2022.csv", col_names = TRUE)

PRIMAPHist %<>%
  pivot_longer(!c(1:6),names_to="Year", values_to = "Value")

#Some additional pre-processing for 0 and NA values
PRIMAPHist <- PRIMAPHist[PRIMAPHist$Value != 0,]

PRIMAPHist %<>%
  na.omit(PRIMAPHist$Value)

#No readable description available. Creating...
PRIMAPHist %<>%
  mutate(Description = case_when(entity == "CH4" ~ "Methane",
                                 entity == "CO2" ~ "Carbon Dioxide",
                                 entity == "N2O" ~ "Nitrous Oxide",
                                 entity == "HFCS (SARGWP100)" ~ "Hydrofluorocarbons (SAR)",
                                 entity == "HFCS (AR4GWP100)" ~ "Hydrofluorocarbons (AR4)",
                                 entity == "PFCS SARGWP100" ~ "Perfluorocarbons (SAR)",
                                 entity == "PFCS (SARGWP100)" ~ "Perfluorocarbons (SAR)",
                                 entity == "PFCS (AR4GWP100)" ~ "Perfluorocarbons (AR4)",
                                 entity == "SF6" ~ "Sulfur Hexafluoride",
                                 entity == "NF3" ~ "Nitrogen Trifluoride",
                                 # entity == "FGASES SARGWP100" ~ "Fluorinated Gases (SAR): HFCs, PFCs, SF$_6$, NF$_3$",
                                 # entity == "FGASES (AR4GWP100)" ~ "Fluorinated Gases (AR4): HFCs, PFCs, SF$_6$, NF$_3$",
                                 entity == "FGASES SARGWP100" ~ "Fluorinated Gases (SAR)",
                                 entity == "FGASES (SARGWP100)" ~ "Fluorinated Gases (SAR)",
                                 entity == "FGASES (AR4GWP100)" ~ "Fluorinated Gases (AR4)",
                                 entity == "KYOTOGHG SARGWP100" ~ "Kyoto greenhouse gases (SAR)",
                                 entity == "KYOTOGHG (SARGWP100)" ~ "Kyoto greenhouse gases (SAR)",
                                 entity == "KYOTOGHG (AR4GWP100)" ~ "Kyoto greenhouse gases (AR4)"),
         SectorDescription = case_when(`category (IPCC2006_PRIMAP)` == "M.0.EL" ~ "National Total excluding LULUCF",
                                       `category (IPCC2006_PRIMAP)` == "1" ~ "Energy",
                                       `category (IPCC2006_PRIMAP)` == "1.A" ~ "Energy - Fuel Combustion Activities",
                                       `category (IPCC2006_PRIMAP)` == "1.B" ~ "Energy - Fugitive Emissions from Fuels",
                                       `category (IPCC2006_PRIMAP)` == "1.B.1" ~ "Energy - Fugitive Emissions - Solid Fuels",
                                       `category (IPCC2006_PRIMAP)` == "1.B.2" ~ "Energy - Fugitive Emissions - Oil and Natural Gas",
                                       `category (IPCC2006_PRIMAP)` == "1.B.3" ~ "Energy - Other Emissions from Energy Production",
                                       `category (IPCC2006_PRIMAP)` == "1.C" ~ "Energy - Carbon Dioxide Transport and Storage",
                                       `category (IPCC2006_PRIMAP)` == "2" ~ "Industrial Processes and Product Use (IPPU)",
                                       `category (IPCC2006_PRIMAP)` == "2.A" ~ "Mineral Industry",
                                       `category (IPCC2006_PRIMAP)` == "2.B" ~ "Chemical Industry",
                                       `category (IPCC2006_PRIMAP)` == "2.C" ~ "Metal Industry",
                                       `category (IPCC2006_PRIMAP)` == "2.D" ~ "Non-Energy Products from Fuels and Solvent Use",
                                       `category (IPCC2006_PRIMAP)` == "2.E" ~ "Electronics Industry",
                                       `category (IPCC2006_PRIMAP)` == "2 F" ~ "Product uses as Substitutes for Ozone Depleting Substances",
                                       `category (IPCC2006_PRIMAP)` == "2.G" ~ "Other Product Manufacture and Use",
                                       `category (IPCC2006_PRIMAP)` == "2.H" ~ "Other Industrial Processes and Product Use",
                                       `category (IPCC2006_PRIMAP)` == "M.AG" ~ "Agriculture, sum of 3.A and M.AG.ELV",
                                       `category (IPCC2006_PRIMAP)` == "3.A" ~ "Livestock",
                                       `category (IPCC2006_PRIMAP)` == "M.AG.ELV" ~ "Agriculture excluding Livestock",
                                       `category (IPCC2006_PRIMAP)` == "4" ~ "Waste",
                                       `category (IPCC2006_PRIMAP)` == "5" ~ "Other"))

#Appears there are 2 different measures for the same type of gas. Excluding SAR measure for consistency. Kyoto greenhouse gases is also an aggregation so excluding to avoid double up (KYOTOGHG (AR4GWP100)).
ExcludeEntityList <- c("HFCS (SARGWP100)", "PFCS SARGWP100", "FGASES SARGWP100", "KYOTOGHG SARGWP100", "KYOTOGHG (AR4GWP100)", "KYOTOGHG (SARGWP100)", "FGASES (SARGWP100)", "PFCS (SARGWP100)")

#Filter by Global figures, exclude above list, use aggregate (National Total excluding LULUCF aka M.0.EL) and use the Country Reported data (HISTCR vs HISTTP "benefit of the doubt")
PRIMAPHist_Earth <- PRIMAPHist[PRIMAPHist$`area (ISO3)` == "EARTH" & PRIMAPHist$`Year` >= 1950 & !(PRIMAPHist$entity %in% ExcludeEntityList) & PRIMAPHist$`category (IPCC2006_PRIMAP)` == "M.0.EL" & PRIMAPHist$`scenario (PRIMAP-hist)` == "HISTCR",]

PRIMAPHist_Earth %<>%
  mutate(Date = as.Date(paste("12/31/", Year), format = "%m/%d/%Y"),
         Description = as.factor(Description))


PRIMAPHist_p1 <- PRIMAPHist_Earth %>%
  plot_ly(
    y = ~Value*1000, 
    x = ~Year, 
    color = ~Description, 
    # text = ~paste(round(Value/sum(Value)*100,2),"%"),
    # text = ~Value,
    text = paste("<b>Emission Type</b>: ", PRIMAPHist_Earth$Description,
                 "<br><b>Year</b>: ", PRIMAPHist_Earth$Year,
                 "<br><b>Tons of CO2</b>: ", round(PRIMAPHist_Earth$Value*1000,2)),
    hoverinfo = "text",
    type = 'scatter', 
    mode = 'lines',
    transforms = list(
    list(
      type = 'aggregate',
      groups = ~Year,
      aggregations = list(
        list(
          target = 'y', func = 'sum', enabled = T
        )
      )
    )
  )
  ) %>%
  ungroup()
  
PRIMAPHist_p1 <- PRIMAPHist_p1 %>% 
  layout(title =list(text = "<b>Global carbon dioxide emissions dominance over time (1950 to 2021)</b>", font = list(size = 19), x = 0.01, yanchor =  'top'),
         yaxis = list(zeroline = FALSE, showgrid = TRUE, nticks = 10, title = "<b>Tonnes</b>"),
         xaxis = list(zeroline = FALSE, showgrid = FALSE, nticks = 15, title = "<b>Year</b>"),
         showlegend = TRUE,
  # legend = list(x = 0.05, y = 0.95))
  legend = list(title = list(text = "<b> Emission type </b>"),
                orientation = "h",   # show entries horizontally
                xanchor = "center",  # use center of legend as anchor
                x = 0.5,             # put legend in center of x-axis
                y = -0.2))

PRIMAPHist_p1

#Too many colours to help with colour blindness but each bar is labelled through hover over anyway with no other dimensions added. May be worth removing legend? But left in to draw immediate attention to Carbon Dioxide.

#Unfortunately had issues with annotations... so could not directly label without causing weird skewing?

```
> **Data Source**: Gütschow, Johannes & Pflüger, Mika 2022, *“The PRIMAP-hist national historical emissions time series (1750-2021) v2.4”*. Zenodo website, https://zenodo.org/record/7179775.<br>
> **SAR vs AR4 use reference**: Gillenwater, M 2010, updated 2022, *"What is a Global warming potential? And which one do I use?"*. GHG Management Institute website, https://ghginstitute.org/2010/06/28/what-is-a-global-warming-potential/<br>

***

* There is always a limit to how much of these that can be absorbed, and unfortunately year on year we can see Carbon Dioxide emissions grow.
<br><br>
So the next question is - now that we've identified the dominating emission, who are the largest emitters of carbon dioxide?
<br>And where does Australia sit in comparison to other countries?

### Who are the largest emitters by outright amount? <br>(5/8)

```{r fig.width = 21, fig.height = 7, fig.align = "center", echo=FALSE}

GlobalMTCO2 <- read_csv("GCB2021v34_MtCO2_flat.csv", col_names = TRUE)

GlobalMTCO2 %<>%
  pivot_longer(!c(1:3),names_to="Type", values_to = "MtCO2")

#Some additional pre-processing for 0 and NA values
GlobalMTCO2 <- GlobalMTCO2[GlobalMTCO2$MtCO2 != 0,]

GlobalMTCO2 %<>%
  na.omit(GlobalMTCO2$MtCO2)

#Some issues with the country names, have to re-source them based on ISO3 code.
GlobalMTCO2$CountryName <- countrycode(GlobalMTCO2$`ISO 3166-1 alpha-3`, "iso3c", "iso.name.en")
GlobalMTCO2$Continent <- countrycode(GlobalMTCO2$`ISO 3166-1 alpha-3`, "iso3c", "continent")

AmountTop20 <- GlobalMTCO2[GlobalMTCO2$`ISO 3166-1 alpha-3` != "WLD" & GlobalMTCO2$`Year` == 2020 & GlobalMTCO2$`Type` == "Total",]
AmountTop20 %<>% arrange(-MtCO2)
AmountTop20 <- AmountTop20[1:20,]

MTCO2_p2 <- plot_ly(data = AmountTop20,
                    y = ~CountryName, 
                    x = ~MtCO2*1000000, 
                    color = ~Continent,
                    colors = c("#a6cee3","#1f78b4","#b2df8a","#33a02c","#fb9a99"), #Using https://colorbrewer2.org/?type=qualitative&scheme=Paired&n=5, potential colour blindness issues but best effort given the number of qualitative values
                    # frame = ~Year, 
                    text = paste("<b>Country</b>: ", AmountTop20$CountryName,
                                 "<br><b>Continent</b>: ", AmountTop20$Continent,
                                 "<br><b>Year</b>: ", AmountTop20$Year,
                                 "<br><b>Tons of CO2</b>: ", round(AmountTop20$MtCO2*1000000,2)),
                    # hoverinfo = "text",
                    textposition = 'none',
                    type = 'bar')
  
MTCO2_p2 <- MTCO2_p2 %>% layout(margin = list(l = -10, r = 0, t = 50, b = 0, pad = 5),
  title =list(text = "<b>China dominates carbon dioxide emissions in 2020</b>", font = list(size = 19), x = 0.01, yanchor =  'top'),
                                yaxis = list(zeroline = FALSE, showgrid = FALSE, title = "", categoryorder = "total ascending"),
                                xaxis = list(zeroline = FALSE, showgrid = TRUE, nticks = 15, title = "<b>Tons of Carbon Dioxide Emitted</b>"),
                                showlegend = TRUE,
  # legend = list(x = 0.8, y = 0.5))
  legend = list(title = list(text = "<b>Continent</b>"),
                orientation = "h",   # show entries horizontally
                xanchor = "center",  # use center of legend as anchor
                x = 0.5,             # put legend in center of x-axis
                y = -0.2))

MTCO2_p2 %>%
  # animation_opts(1000, redraw = TRUE) %>%
  # animation_slider(currentvalue = list(prefix = "Year: ", font = list(color="black"))) %>% 
  add_annotations(x=500*1000000,
                  y="Australia",
                  xref = "x",
                  yref = "y",
                  text = "<i><b>Australia sits at Rank #16</b></i>",
                  xanchor = 'left',
                  showarrow = T,
                  ax = 300,
                  ay = 0)

# Issues with annotations not updating correctly and interfering with colour
# Is it misleading to only take a subset of countries? Removed animation as a result.

```
> **Data Source**: Andrew, Robbie M. & Peters, Glen P. 2021, *“The Global Carbon Project's fossil CO2 emissions dataset”*, Zenodo website, https://zenodo.org/record/5569235  <br>
> **Colour scheme reference:**: Brewer c, Harrower M, and The Pennsylvania State University (2021). *ColorBrewer: Color Advice for Maps.*. Accessed 3 December, 2022, Published by Color Brewer 2 website: https://colorbrewer2.org/?type=qualitative&scheme=Paired&n=5

***

* But Australia is such a large country with a lot of natural resources and comparatively - a much smaller population than some of these countries. So it only makes sense that Australia would be here...
* What if we take emissions per person into account?

### What about when population is taken into consideration? <br>(6/8)

```{r fig.width = 21, fig.height = 7, fig.align = "center", echo=FALSE}

PerCapitaTop20 <- GlobalMTCO2[GlobalMTCO2$`ISO 3166-1 alpha-3` != "WLD" & GlobalMTCO2$`Year` == 2020 & GlobalMTCO2$`Type` == "Per Capita",]
PerCapitaTop20 %<>% arrange(-MtCO2)
PerCapitaTop20 <- PerCapitaTop20[1:20,]

MTCO2_p3 <- plot_ly(data = PerCapitaTop20,
    y = ~CountryName, 
    x = ~MtCO2, 
    color = ~Continent,
    colors = c("#a6cee3","#1f78b4","#b2df8a","#33a02c"), #Using https://colorbrewer2.org/?type=qualitative&scheme=Paired&n=4, potential colour blindness issues but best effort given the number of qualitative values
    # frame = ~Year, 
    text = paste("<b>Country</b>: ", PerCapitaTop20$CountryName,
                 "<br><b>Continent</b>: ", PerCapitaTop20$Continent,
                 "<br><b>Year</b>: ", PerCapitaTop20$Year,
                 "<br><b>Per Capita Tons of CO2</b>: ", round(PerCapitaTop20$MtCO2,2)),
    # hoverinfo = "text",
    type = 'bar')
  
MTCO2_p3 <- MTCO2_p3 %>% layout(margin = list(l = -10, r = 0, t = 50, b = 0, pad = 5),
                                title =list(text = "<b>Qatar leads carbon dioxide emissions per capita in 2020</b>", font = list(size = 19), x = 0.01, yanchor =  'top'),
                                yaxis = list(zeroline = FALSE, showgrid = FALSE, title = "", categoryorder = "total ascending"),
                                xaxis = list(zeroline = FALSE, showgrid = TRUE, nticks = 15, title = "<b>Per capita tons of carbon dioxide Emitted</b>"),
                                showlegend = TRUE,
  legend = list(title = list(text = "<b> Continent </b>"),
                orientation = "h",   # show entries horizontally
                xanchor = "center",  # use center of legend as anchor
                x = 0.5,             # put legend in center of x-axis
                y = -0.2))

MTCO2_p3 %>%
  # animation_opts(1000, redraw = TRUE) %>%
  # animation_slider(currentvalue = list(prefix = "Year: ", font = list(color="black"))) %>% 
  add_annotations(x=15.5,
                  y="Australia",
                  xref = "x",
                  yref = "y",
                  text = "<i><b>Australia sits at Rank #11</i></b>",
                  xanchor = 'left',
                  showarrow = T,
                  ax = 80,
                  ay = 0)

#Removing animation as it could be deceptive to only see a subset of countries over the years

```
> **Data Source**: Andrew, Robbie M. & Peters, Glen P. 2021, *“The Global Carbon Project's fossil CO2 emissions dataset”*, Zenodo website, https://zenodo.org/record/5569235  <br>
> **Colour scheme reference:**: Brewer c, Harrower M, and The Pennsylvania State University (2021). *ColorBrewer: Color Advice for Maps.*. Accessed 3 December, 2022, Published by Color Brewer 2 website: https://colorbrewer2.org/?type=qualitative&scheme=Paired&n=4

***

* Unfortunately Australia ranks even *higher* with emissions per person taken into view.
* So as we're diving deeper into this, we need to find out the root cause of these emissions in Australia.

### Where is all of Australia's Carbon Dioxide coming from? <br>(7/8)

```{r fig.width = 21, fig.height = 7, fig.align = "center", echo=FALSE}

#Exclude aggregated sectors, focus on Australia in 2020 for CO2 only.
PRIMAPHist_Aus <- PRIMAPHist[PRIMAPHist$`area (ISO3)` == "AUS" & PRIMAPHist$`Year` == 2020 & !(PRIMAPHist$entity %in% ExcludeEntityList) & PRIMAPHist$`scenario (PRIMAP-hist)` == "HISTCR" & PRIMAPHist$`category (IPCC2006_PRIMAP)` != "M.0.EL" & PRIMAPHist$`category (IPCC2006_PRIMAP)` != "M.AG" & PRIMAPHist$`category (IPCC2006_PRIMAP)` != "1" & PRIMAPHist$`category (IPCC2006_PRIMAP)` != "1.B" & PRIMAPHist$`category (IPCC2006_PRIMAP)` != "2" & PRIMAPHist$entity == "CO2",]

PRIMAPHist_Aus %<>%
  mutate(Date = as.Date(paste("12/31/", Year), format = "%m/%d/%Y"),
         Description = as.factor(Description),
         SectorDescription = as.factor(SectorDescription),
         Prop = Value/sum(Value))

PRIMAPHist_p2 <- plot_ly(data = PRIMAPHist_Aus,
                    y = ~SectorDescription, 
                    x = ~Value*1000, 
                    color = ~SectorDescription,
                    # frame = ~Year, 
                    text = paste("<b>Country</b>: Australia",
                                 "<br><b>Category/Sector</b>: ", PRIMAPHist_Aus$SectorDescription,
                                 "<br><b>Year</b>: ", PRIMAPHist_Aus$Year,
                                 "<br><b>Tons</b>: ", round(PRIMAPHist_Aus$Value*1000,2),
                                 "<br><b>Percentage of total (%)</b>: ", round(PRIMAPHist_Aus$Prop*100,2)),
                    # hoverinfo = "text",
                    textposition = 'none',
                    type = 'bar')
  
PRIMAPHist_p2 <- PRIMAPHist_p2 %>% layout(margin = list(l = -10, r = 0, t = 70, b = -5, pad = 5),
                                          title =list(text = "<b>Energy production via fuel combustion is the clear segment responsible \nfor carbon dioxide emissions in Australia for 2020</b>", font = list(size = 19), x = 0.01, yanchor =  'top'),
                                          yaxis = list(zeroline = FALSE, showgrid = FALSE, title = "<b>Category/Sector</b>", categoryorder = "total ascending"),
                                          xaxis = list(zeroline = FALSE, showgrid = TRUE, nticks = 15, title = "<b>Tons of Carbon Dioxide Emissions</b>"),
                                          showlegend = FALSE,
                                          barmode = "stack")

PRIMAPHist_p2 
# %>%
  # animation_opts(1000, redraw = TRUE) %>%
  # animation_slider(currentvalue = list(prefix = "Year: ", font = list(color="black"))) %>% 
  # add_annotations(x=500,
  #                 y="Australia",
  #                 xref = "x",
  #                 yref = "y",
  #                 text = "Australia sits at Rank #16",
  #                 xanchor = 'left',
  #                 showarrow = T,
  #                 ax = 300,
  #                 ay = 0)

# Issues with annotations not updating correctly and interfering with colour
# Is it misleading to only take a subset of countries? Removed animation as a result.

# write.csv(as.data.frame(PRIMAPHist_Aus), "PRIMAPHIST - AUS Review Aggregates.csv", row.names=T)

#Should this be cumulative? Will 1 year in isolation be ignoring "seasonal" type trends?
#Leaving as is, as cumulative figures might "colour" any recent efforts to improve emissions.

#Too many colours to help with colour blindness but each bar is labelled anyway with no other dimensions added.

```
> **Data Source**: Gütschow, Johannes & Pflüger, Mika 2022, *“The PRIMAP-hist national historical emissions time series (1750-2021) v2.4”*. Zenodo website, https://zenodo.org/record/7179775.<br><br>

***

**Important Note**: Using the displayed data source reveals these findings, however results from the National Greenhouse Gas Inventory (used by the CSIRO) still shows Energy as a leading source - the breakdown shows an entirely different make up and categorisation. Would recommend comparing other sources for more insight. Suspect there may be some level of estimation involved regardless of the source referenced. Refer to link below for more information. <br><br>

**National Greenhouse Gas reference**: CSIRO (Australia’s National Science Agency), *‘What are the sources of Australia's greenhouse gases?’*. Accessed 3 December 2022, CSIRO website: https://www.csiro.au/en/research/environmental-impacts/climate-change/climate-change-qa/sources-of-ghg-gases


### If Energy is the biggest problem - what have we been doing to turn this around? <br>(8/8)

```{r fig.width = 21, fig.height = 7, fig.align = "center", echo=FALSE}

# AES2021Aus <- read.xlsx("Australian Energy Statistics 2021 Table C_0.xlsx", sheet="AUS", colNames = TRUE, startRow = 5)
# 
# AES2021Aus %<>%  
#   filter(!row_number() %in% c(1, 62, 63, 64))
# 
# colnames(AES2021Aus)[1] <- 'FY'
# 
# AES2021Aus %<>%
#   mutate(across(where(is.character), ~na_if(., "na")),
#          Coal = as.numeric(Coal),
#          Oil = as.numeric(Oil),
#          Gas = as.numeric(Gas),
#          Renewables = as.numeric(Renewables),
#          Statistical.Discrepancy.b = as.numeric(Statistical.Discrepancy.b),
#          Total = as.numeric(Total),
#          RenewableProp = Renewables/(Coal+Oil+Gas+Renewables))
#Realised that production may outweight consumption...

AES2021Aus <- read.xlsx("Australian Energy Statistics 2021 Table O_0.xlsx", sheet="AUS", colNames = TRUE, startRow = 5)

AES2021Aus %<>%
  filter(!row_number() %in% c(1, 2, 9, 19:38))

colnames(AES2021Aus)[1] <- 'Fuel Type'

AES2021Aus %<>%
  pivot_longer(!c(1),names_to="Financial Year", values_to = "GWh")

AES2021Aus %<>%
  mutate(GWh = as.numeric(GWh))

AES2021Aus %<>%
  pivot_wider(names_from = `Fuel Type`, values_from = GWh)

#Are the other dataframes in a tidy format like this? May need to check and rework...

AES2021Aus %<>%
  mutate(RenewablePortion = `Total renewable`/(`Total renewable`+`Total non-renewable`),
         NonRenewablePortion = `Total non-renewable`/(`Total renewable`+`Total non-renewable`))

# AES2021Aus_p1 <- plot_ly(data = AES2021Aus, x = ~`Financial Year`, y = ~`Total renewable`, type = 'bar', name = 'Renewable energy (GWh)') %>% 
#   add_trace(y = ~`Total non-renewable`, name = 'Non-Renewable energy (GWh)') %>% 
#   layout(yaxis = list(title = 'GWh'), barmode = 'stack')

AES2021Aus_p1 <- plot_ly(data = AES2021Aus, 
                         x = ~`Financial Year`, 
                         y = ~RenewablePortion*100,
                         marker = list(color = "#1b9e77"),
                         type = 'bar', 
                         name = 'Renewable energy (% of total)',
                         text = paste("<b>Financial Year</b>: ", AES2021Aus$`Financial Year`,
                                      "<br><b>Renewable energy portion</b>: ", round(AES2021Aus$RenewablePortion*100,2),"%",
                                      "<br><b>Renewable energy (GWh)</b>: ", round(AES2021Aus$`Total renewable`,2)),
                         textposition = 'none') %>%
  add_trace(y = ~NonRenewablePortion*100, 
            marker = list(color = "#7570b3"),
            name = 'Non-Renewable energy (% of total)',
            text = paste("<b>Financial Year</b>: ", AES2021Aus$`Financial Year`,
                         "<br><b>Non-Renewable energy portion</b>: ", round(AES2021Aus$NonRenewablePortion*100,2),"%",
                         "<br><b>Non-Renewable energy (GWh)</b>: ", round(AES2021Aus$`Total non-renewable`,2)),
                         textposition = 'none') %>%
  layout(margin = list(l = -10, r = 20, t = 50, b = -5, pad = 5),
         title =list(text = "<b>Australia demonstrates a slow transition to renewable energy generation</b>", font = list(size = 19), x = 0.01, yanchor =  'top'),
         yaxis = list(zeroline = FALSE, showgrid = FALSE, title = "<b>Percentage of total energy generated (%)</b>"),
         xaxis = list(zeroline = FALSE, showgrid = TRUE, nticks = 12, title = "<b>Financial year (FY)</b><br>Note: Data is available for calendar year but will only be estimates."),
         barmode = 'stack',
         # legend = list(y = 100, x = 0.65))
         legend = list(title = list(text = "<b>Energy type</b>"),
                orientation = "h",   # show entries horizontally
                xanchor = "center",  # use center of legend as anchor
                x = 0.5,             # put legend in center of x-axis
                y = -0.3))


AES2021Aus_p1

#Need to note estimates where possible from original data source.

```
> **Data Source**: Australian Government, Department of Climate Change, Energy, the Environment and Water 2021, *“TAustralian Energy Update 2021”*. Department of Climate Change, Energy, the Environment and Water website, https://www.energy.gov.au/publications/australian-energy-update-2021. <br>
> **Colour scheme reference:**: Brewer c, Harrower M, and The Pennsylvania State University (2021). *ColorBrewer: Color Advice for Maps.*. Accessed 3 December, 2022, Published by Color Brewer 2 website: https://colorbrewer2.org/#type=qualitative&scheme=Dark2&n=3

***

* Another factor not displayed, are energy demands year on year. However from a brief review, the amount of GWh generated from non-renewables appears to be in slow decline.
* The positive note is that there is a shift to renewables - but whether this is too late is another thing.

3. Conclusion
================================================================

Column {data-width="1%"}
-----------------------------------------------------------------------

Column {data-width="25%"}
-----------------------------------------------------------------------

```{r fig.cap = "Source: Eildon Power Station - by Benjamin Johnson", echo=FALSE} 

knitr::include_graphics("_DSC3997.jpg")

```

Column {data-width="3%"}
-----------------------------------------------------------------------

Column {data-width="70%"}
-----------------------------------------------------------------------
<br>

To recap, we've found...

1. Global temperatures have been increasing, with the highest global average being **14.78** degrees celsius (with varied increases depending on which country we look at). <br>
Globally this is **0.81 degrees higher than the 1951-1980 average**.
2. **Carbon Dioxide emissions sat at 36.3 billion tonnes for 2021** with the next closest emission only sitting at 1.1 billion tonnes. The rate of emissions are only increasing.
3. Of all the countries for outright carbon dioxide amount - **China emits the most (#1), sitting at 10.6 billion tons for 2020**. Australia ranks #16 at 391 million tons.
4. For carbon dioxide emissions per capita also **in 2020, Qatar leads at 37 per capita tons**. Australia ranks #11 at just over 15 per capita tons.
5. In Australia, 90.7% of carbon dioxide emissions is attributable to Energy - Fuel Combustion activities (Open query. See *Note* for CSIRO view).
6. **In FY20, Australia is only generating 22.6% of it's energy from renewable sources**. But if trends continue, hopefully increases over time.

*Note*: According to the CSIRO [National Greenhouse Gas Inventory 2020] while Energy is still the largest contributor to Carbon Dioxide emissions in Australia, the make up is instead as follows - 

* Energy (burning fossil fuels to produce electricity) contributed 33.6 per cent of the total emissions
* Stationary energy (including manufacturing, mining, residential and commercial fuel use) 20.4 per cent
* Transport 17.6 per cent
* Agriculture 14.6 per cent
* Fugitive emissions 10.0 per cent
* Industrial processes 6.2 per cent
* Waste 2.7 per cent.

<br>
To close off:

* It is important to note, that while we've focused on the biggest contributors - there are also other factors as well to consider as we've found with the differing categorisations.
* We have established that our reliance on fossil fuels for energy production is a significant contributor to global warming.
* Australia's transition to renewable energy has a long way to go but progress is being made.
* From here, there is a clear need to accelerate this transition to support the global effort to slow down climate change as much as possible so Australia does not see continuing increases from the **1.47** increase which is already our reality today.
<br>

***

**National Greenhouse Gas reference**: CSIRO (Australia’s National Science Agency), *‘What are the sources of Australia's greenhouse gases?’*. Accessed 3 December 2022, CSIRO website: https://www.csiro.au/en/research/environmental-impacts/climate-change/climate-change-qa/sources-of-ghg-gases

Column {data-width="1%"}
-----------------------------------------------------------------------

```{r fig.width = 21, fig.height = 21, fig.align = "center", echo=FALSE}
# Appendix - Obsolete code

#Point map of global cities, not fit for purpose when displaying country temperatures...
# global_cities <- read_csv("worldcities.csv", col_names = TRUE)
# 
# 
# map <- global_cities %>% 
#   sample_n(50) %>% 
#   leaflet() %>% 
#   addTiles() %>%
#   addMarkers(~lng, ~lat, popup = ~htmlEscape(paste(city,":", iso3 )))
# 
# map

#Warning: Data and computation intensive. Can use CSV to skip processing. 
#Ended up using the extract from Antti Lipponen due to complexity with mapping lat long co-ordinates and aggregating to country boundaries...

# CCKP_MeanTemp <- nc_open("timeseries-tas-annual-mean_cru_annual_cru-ts4.06-timeseries_mean_1901-2021.nc")
# print(CCKP_MeanTemp)
# 
# lon <- ncvar_get(CCKP_MeanTemp, "lon")
# lat <- ncvar_get(CCKP_MeanTemp, "lat")
# time <- ncvar_get(CCKP_MeanTemp, "time")
# 
# head(time)
# tunits <- ncatt_get(CCKP_MeanTemp, "time", "units")
# 
# degC_array <- ncvar_get(CCKP_MeanTemp, "timeseries-tas-annual-mean")
# fillvalue <- ncatt_get(CCKP_MeanTemp, "timeseries-tas-annual-mean", "_FillValue")
# 
# degC_array[degC_array==fillvalue$value] <- NA
# #degC_array
# 
# print(tunits)
# 
# tustr <- strsplit(tunits$value, " ")
# tdstr <- strsplit(unlist(tustr)[3], "-")
# tmonth <- as.integer(unlist(tdstr)[2])
# tday <- as.integer(unlist(tdstr)[3])
# tyear <- as.integer(unlist(tdstr)[1])
# time_obs <- as.POSIXct(chron(time,origin=c(tmonth, tday, tyear)))
# 
# #degC_slice <- degC_array[ , , 25]
# #image(lon, lat, degC_slice)
# 
# lonlattime <- as.matrix(expand.grid(lon,lat,time_obs))
# degC_vec_long <- as.vector(degC_array)
# length(degC_vec_long)
# 
# degC_obs <- data.frame(cbind(lonlattime,degC_vec_long))
# #head(degC_obs)
# 
# dim(degC_obs)
# 
# CCKP_degC <- na.omit(degC_obs)
# colnames(CCKP_degC) <- c("Long","Lat","Date","Mean_Temp")
# 
# #write.csv(as.data.frame(CCKP_degC), "timeseries-tas-annual-mean_cru_annual_cru-ts4.06-timeseries_mean_1901-2021.csv", row.names=T)
# 
# plot_mapbox(data = CCKP_degC[CCKP_degC$Date == "2021-01-16 11:00:00",]) %>%
#   add_markers(
#     x = ~Long, 
#     y = ~Lat, 
#     color = ~Mean_Temp,
#     colors = "Accent",
#     text = ~Mean_Temp,
#     hoverinfo = "text"
#   )

#> Reference: https://towardsdatascience.com/how-to-crack-open-netcdf-files-in-r-and-extract-data-as-time-series-24107b70dcd
#> Reference: https://rpubs.com/boyerag/297592

#Removed as it does not add anything to the overarching story beyond next level detail

### Continent view
# install.packages("geodata")
# library(geodata)

# headers = read.csv("tas_timeseries_annual_cru_1901-2021_AUS.csv", skip = 1, header = F, nrows = 1, as.is = T)
# AnnualAUTemps = read.csv("tas_timeseries_annual_cru_1901-2021_AUS.csv", skip = 2, header = F)
# colnames(AnnualAUTemps)= headers
# colnames(AnnualAUTemps)[1] <- "Year"
# 
# AnnualAUTemps %<>%
#   pivot_longer(!c(1),names_to="Region/State", values_to = "Temp")
# 
# ausData <- getData("GADM", country = "AUS", level = 1, download = TRUE) #Get state level data
# 
# aus_map <- leaflet(ausData, options=leafletOptions(worldCopyJump = TRUE, minZoom = 1.5, zoomSnap = 0.25)) %>%
#   setView(134, -29, zoom=4)
# 
# aus_map %>%
#   addPolygons(
#   # fillColor = ~pal(Temp),
#   popup = paste("State: ",ausData$NAME_1),
#   weight = 0.7,
#   opacity = 1,
#   color = "black",
#   dashArray = "3",
#   fillOpacity = 0.5,
#   smoothFactor = 0.2,
#   stroke = TRUE,
#   highlightOptions = highlightOptions(
#     weight = 3,
#     color = "#666",
#     dashArray = "",
#     fillOpacity = 0.5,
#     bringToFront = TRUE),
#   # label = labels,
#   labelOptions = labelOptions(
#     style = list("font-weight" = "normal", padding = "3px 8px"),
#     textsize = "15px",
#     direction = "auto"),
#   group = "AU")
# %>%
#   addLegend(pal = pal,
#             values = ~Temp,
#             opacity = 0.7,
#             title = "Mean Temperature Anomaly",
#             labFormat = labelFormat(suffix = "°C"),
#             position = "bottomleft")

#> Reference for Aus data: https://www.youtube.com/watch?v=AjFVQ8ujoiE

#Tried to include valueBox in the commentary section unsuccessfully. Strips it back to text.
#flexdashboard::valueBox(value = paste("+",GlobalTempAnomaly[GlobalTempAnomaly$end_year == 2019,]$mean_temp_anomaly), caption = "Temperature difference from 1951-1980 average", icon = "fa-comments")

#Unused references but useful for further research/insight

# > BP source: https://www.bp.com/en/global/corporate/energy-economics/statistical-review-of-world-energy/primary-energy.html<br>
# > Ember Electricity Source: https://ember-climate.org/data-catalogue/yearly-electricity-data/<br>
# > OurWorldinData Source: https://ourworldindata.org/electricity-mix<br>
# > Australian Energy Update 2021: https://www.energy.gov.au/publications/australian-energy-update-2021<br>
# > Australian Energy by Fuel Type: https://www.energy.gov.au/publications/australian-energy-statistics-table-o-electricity-generation-fuel-type-2020-21-and-2021<br>
# > Energy Data links: https://www.energy.gov.au/government-priorities/energy-data<br>
# > World Bank Group - Climate Change Knowledge Portal (World Bank Group, Climate Change Knowledge Portal.): https://climateknowledgeportal.worldbank.org/download-data
# > GeoJSON Source: https://geojson-maps.ash.ms/<br>
# > World Cities: https://simplemaps.com/data/world-cities<br>
# > Leaflet reference: https://rstudio.github.io/leaflet/json.html<br>
# > For Zoom adjustmnet reference: https://leafletjs.com/examples/zoom-levels/<br>
# > Alt Reference: https://globalecoguy.org/the-three-most-important-graphs-in-climate-change-e64d3f4ed76<br>
# > Potential data source: https://www.climatewatchdata.org/ghg-emissions?breakBy=gas&chartType=area&end_year=2019&regions=WORLD&source=CAIT&start_year=1990<br>
# > Countrycode reference: https://cran.r-project.org/web/packages/countrycode/countrycode.pdf<br>
# > Plotly reference: https://plotly.com/r/reference/#layout-annotations<br>

```