```{r, include = FALSE}
# Load libraries 
library(tidyr)
library(dplyr)
library(readr)
library(ggplot2)
library(lubridate)
library(plotly)
library(sf)
library(crosstalk)
library(cartogram)
library(ggthemes)
```

## About the data

[Indonesia earthquakes datasets](https://github.com/horankev/quake_data)

These datasets contain spatio-temporal information about earthquakes in Indonesia.

`quakes_df.rds` contains records for individual earthquakes. `provinces_df.rds` contains records for each province, including aggregations of earthquake counts and densities.

## Steps to take

::: {style="font-size: 80%;"}

1. Read in the two data sets
2. Subset the data to contain only quakes in November 2023, small enough to experiment with
3. Make a map with earthquake locations plotted
4. Make an interactive plot of magnitude and date
5. Link the two plots so that a user can select points in one plot and examine where they are in the other
6. Develop an animation of earthquakes over time
7. Create a choropleth map of earthquakes by province
8. Convert the map to a cartogram
9. Add interactive mouse-over
10. Create a quarto document to start a report on earthquakes in Indonesia

:::

## 1. Read

```{r}
#| eval: false
#| echo: true
quakes_sf <- readRDS("data/earthquakes/quakes_df.rds")
provinces_sf <- readRDS("data/earthquakes/provinces_df.rds")

# Check map
ggplot(provinces_sf) + geom_sf()
```

## 2. Subset the data

```{r}
#| eval: false
#| echo: true
quakes_df <- st_drop_geometry(quakes_sf)
quakes_df <- quakes_df |>
  bind_cols(st_coordinates(quakes_sf)) 
quakes_11_2023 <- quakes_df |>
  filter(year(ymd) == 2023, month(ymd) == 11)
```

## 3. Make a map

```{r}
#| eval: false
#| echo: true
ggplot(provinces_sf) + 
  geom_sf() +
  geom_point(data=quakes_11_2023, aes(x=X, y=Y)) +
  theme_map()
```

## 4. Make an interactive plot

```{r}
#| eval: false
#| echo: true
p1 <- ggplot(quakes_11_2023, aes(x=ymd, y=mag)) +
  geom_point()
ggplotly(p1)
```

Note: `sf` objects don't work with `plotly` if there are multiple geometry types in them.

```{r}
#| eval: false
#| echo: true
quakes_11_2023_shared <- highlight_key(quakes_11_2023)

p1 <- ggplot(quakes_11_2023_shared, aes(x=ymd, y=mag)) +
  geom_point() +
  theme_minimal()

gp1 <- ggplotly(p1, width=400, height=400) |>
  config(displayModeBar = FALSE) |>
  highlight(on = "plotly_selected", 
            off = "plotly_doubleclick") 

# multiple types of geometry in provinces_sf needs fixing
provinces_sf2 <- st_cast(provinces_sf, "MULTIPOLYGON")
p2 <- ggplot() + 
  geom_sf(data=provinces_sf2) +
  geom_point(data=quakes_11_2023_shared, aes(x=X, y=Y)) +
  theme_map()

gp2 <- ggplotly(p2, width=800, height=400) |>
  config(displayModeBar = FALSE) |>
  highlight(on = "plotly_selected", 
            off = "plotly_doubleclick") 

bscols(gp1, gp2, widths = c(4, 8))
```


## 7. Create a choropleth map

```{r}
#| eval: false
#| echo: true
quakes_2023 <- quakes_df |>
  filter(year(ymd) == 2023) |>
  group_by(province) |>
  summarise(mag = mean(mag, na.rm=TRUE),
            n = n(), na.rm=TRUE)

provinces_sf2 <- left_join(provinces_sf2, quakes_2023)

p3 <- ggplot(provinces_sf2) + 
  geom_sf(aes(fill=mag, label=province)) +
  theme_map() +
  scale_fill_viridis_c()
```
## 8. Convert to a cartogram

Population data can be found [here](https://sultra.bps.go.id/en/statistics-table/1/NDc3OCMx/population-by-province-in-indonesia--thousand---2019--2023.html).

```{r}
#| eval: false
#| echo: true
pop <- read_csv("data/province_pop.csv") |>
  rename(province = Provinsi)
provinces_sf2 <- left_join(provinces_sf2, pop)

indonesia_carto <- provinces_sf2 |> 
  st_transform(3395) |> 
  cartogram_cont("2023") |>
  st_transform("WGS84") 

indonesia_carto <- st_cast(indonesia_carto, "MULTIPOLYGON")

ggplot(indonesia_carto) + geom_sf()

indonesia_carto <- left_join(indonesia_carto, quakes_2023)

p4 <- ggplot(indonesia_carto) + 
  geom_sf(aes(fill=mag, label=province)) +
  theme_map() +
  scale_fill_viridis_c()
```

## 9. Interactivity

```{r}
#| eval: false
#| echo: true

ggplotly(p3)

ggplotly(p4)
```

## 10. Embed your best examples into a quarto document

Once you have your quarto document rendering, change the data to use 2022, and re-render, to see how fast it is to re-do the analysis if the data changes.

## Resources


::: bottom
<span style="display:inline-block;"><a rel="license" href="http://creativecommons.org/licenses/by-nc-sa/4.0/"><img alt="Creative Commons License" style="border-width:0" src="https://i.creativecommons.org/l/by-nc-sa/4.0/88x31.png" /></a></span><span style="display:inline-block;"> This work is licensed under a <a rel="license" href="http://creativecommons.org/licenses/by-nc-sa/4.0/">Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International License</a>.</span>
:::