PlotExercises_report.Rmd

---
title: "Plotting exercises with R"
output:
  html_document:
    highlight: pygments
    code_folding: show
    depth: 4
    number_sections: no
    theme: sandstone
    toc: yes
    toc_float:
      collapsed: yes
      smooth_scroll: yes
date: "`r Sys.Date()`"
---

<style type="text/css">
  body{
  font-family: Exo;
}
</style>

```{r setup, include=FALSE}
knitr::opts_chunk$set(tidy=FALSE,
                      echo = TRUE, 
                      cache = TRUE,
                      message=FALSE,
                      warning=FALSE,
                      fig.align='center')
knitr::opts_knit$set(progress = TRUE, verbose = TRUE)

library(tidyverse)
library(gganimate)
library(here)
library(janitor)
library(hues)
library(ggpubr)
library(magick)
library(sf)
library(giscoR)
library(stars)
library(rayshader)
library(MetBrewer)
library(colorspace)

i_am('PlotExercises_report.Rmd')

# library(showtext)
# font_add_google(name = "Exo", family = "Exo", regular.wt = 400, bold.wt = 700)
# showtext.auto()

library(extrafont)
font_import()
loadfonts()
```

Practice of different types of data visualization 

****

## Bar race animation charts

Following code from https://www.r-bloggers.com/2020/01/how-to-create-bar-race-animation-charts-in-r/.

For this example it will be used a data set from the world data bank for the [ratio of female to male labor force participation rate](https://databank.worldbank.org/source/gender-statistics/Type/TABLE/preview/on#). The data is in the datasets folder.

### Clean the data

```{r}
dt <- read_csv(here('datasets/RatioGenderLaborForce.csv'),show_col_types = FALSE, na = "..") %>%
  select(-"Series Name", -"Series Code" ) %>%
  drop_na() %>%
  mutate_at(vars(contains("YR")),as.numeric) %>%
  pivot_longer(cols = -c("Country Name","Country Code"), names_to = "year") %>%
  janitor::clean_names() %>%
  mutate(year = as.numeric(stringr::str_sub(year,1,4))) 
```

### Format the data in order to keep only the top 10 worse countries for every given year

```{r}
dt_formatted <- dt %>%
  group_by(year) %>% 
  mutate(rank = rank(-value),
         value = round(value,0)) %>%
  group_by(country_name) 

dt_formattedTop <- dt_formatted %>%
  filter(rank <= 10) %>%
  ungroup()

dt_formattedbottom <- dt_formatted %>%
  filter(rank > max(unique(dt_formatted$rank) - 10)) %>%
  ungroup()
```

### Build static plot

```{r}
staticplot1 = ggplot(dt_formattedTop, aes(rank, group = country_name,
                fill = as.factor(country_name), color = as.factor(country_name))) +
  geom_tile(aes(y = value/2,
                height = value,
                width = 0.9), alpha = 0.8, color = NA) +
  geom_text(aes(y = 0, label = paste(country_name, " ")), vjust = 0.2, hjust = 1) +
  geom_text(aes(y=value, label = as.character(round(value,0)), hjust=0)) +
  coord_flip(clip = "off", expand = FALSE) +
  ylim(0,110) +
  scale_x_reverse() +
  scale_fill_iwanthue() +
  scale_color_iwanthue() +
  guides(color = "none", fill = "none") +
theme(axis.line=element_blank(),
      axis.text=element_blank(),
      axis.ticks=element_blank(),
      axis.title=element_blank(),
      legend.position="none",
      panel.background=element_blank(),
      panel.border=element_blank(),
      panel.grid=element_blank(),
      # panel.grid.major.x = element_line( size=.1, color="grey" ),
      # panel.grid.minor.x = element_line( size=.1, color="grey" ),
      plot.subtitle=element_text(hjust=1),
      plot.background=element_blank(), plot.margin = margin(1,1, 1, 4, "cm"),
      text=element_text(family="Exo")) +  
  transition_states(year, transition_length = 3, state_length = 0,
                    wrap = FALSE) +
  view_follow(fixed_x = TRUE, fixed_y = TRUE)  +
  labs(subtitle = "Best 10 Countries : {closest_state}") +
   enter_fade() + 
  exit_fade() +
  ease_aes('linear')

staticplot2 = ggplot(dt_formattedbottom, aes(rank, group = country_name,
                fill = as.factor(country_name), color = as.factor(country_name))) +
  geom_tile(aes(y = value/2,
                height = value,
                width = 0.9), alpha = 0.8, color = NA) +
  geom_text(aes(y = 0, label = paste(country_name, " ")), vjust = 0.2, hjust = 1) +
  geom_text(aes(y=value, label = as.character(round(value,0)), hjust=0)) + #-0.5
  coord_flip(clip = "off", expand = FALSE) +
    ylim(0,110) +
  scale_x_reverse() +
  scale_fill_iwanthue() +
  scale_color_iwanthue() +
  guides(color = "none", fill = "none") +
theme(axis.line=element_blank(),
      axis.text=element_blank(),
      axis.ticks=element_blank(),
      axis.title=element_blank(),
      legend.position="none",
      panel.background=element_blank(),
      panel.border=element_blank(),
      panel.grid=element_blank(),
      # panel.grid.major.x = element_line( size=.1, color="grey" ),
      # panel.grid.minor.x = element_line( size=.1, color="grey" ),
      plot.subtitle=element_text(hjust=1),
      plot.background=element_blank(), plot.margin = margin(1,1, 1, 4, "cm"),
      text=element_text(family="Exo")) +  
  transition_states(year, transition_length = 3, state_length = 0,
                    wrap = FALSE) +
  view_follow(fixed_x = TRUE, fixed_y = TRUE)  +
  labs(subtitle = "Worse 10 Countries : {closest_state}") +
  enter_fade() + 
  exit_fade() +
  ease_aes('linear')

```

### Make and render animation with two plots

```{r}
frames=300
gif1 <- animate(staticplot1, frames, width = 600, height = 350, res = 100, duration = 30,
        renderer = gifski_renderer("outputs/gganim1.gif"), end_pause = 20)
gif2 <- animate(staticplot2, frames,width = 600, height = 350, res = 100, duration = 30,
        renderer = gifski_renderer("outputs/gganim2.gif"), end_pause = 20)

a_mgif <- image_read(gif1)
b_mgif <- image_read(gif2)
if(length(a_mgif) == length(b_mgif)){
  new_gif <- image_append(c(a_mgif[1], b_mgif[1]), stack =TRUE)
  for(i in 2:length(a_mgif)){
    combined <- image_append(c(a_mgif[i], b_mgif[i]), stack =TRUE)
    new_gif <- c(new_gif, combined)
  }
}

gif <- new_gif %>%
  image_annotate("Ratio of female to male labor force participation rate", 
                 size = 16, font = "Exo", gravity = "north") %>%
  image_annotate("Data Source: World Bank Data", 
                 size = 10, font = "Exo", gravity = "southwest")

gif
anim_save("outputs/gganim12.gif",animation = gif)
```

****

## 3D population density density map 

Based on the [visual capitalist](https://www.visualcapitalist.com/cp/population-density-patterns-countries/) post with code from [github](https://github.com/Pecners/kontur_rayshader_tutorial).

Using population data from [kontur](https://www.kontur.io/portfolio/population-dataset/) and country shape data from giscoR to plot the map with the package Rayshader.

[Another post to look](https://spencerschien.info/post/data_viz_how_to/high_quality_rayshader_visuals/) and [this youtube video from the code author](https://www.youtube.com/watch?v=zgFXVhmKNbU).

### Load the data and the shapefile

```{r, eval = FALSE}
data <- st_read(here('datasets/kontur_population_PT_20220630.gpkg'))

#EPSG: 4326 uses a coordinate system on the surface of a sphere or ellipsoid of reference. Think of it as this way: EPSG 4326 uses a coordinate system the same as a GLOBE (curved surface). EPSG 3857 uses a coordinate system the same as a MAP (flat surface)
#"4326": WGS84

pt <- giscoR::gisco_get_countries(
  epsg = "3857", #pseudo-Mercator
  resolution = "3",
  country = "Portugal")  %>% 
  st_transform(crs = st_crs(data))
```

### Limit the data to the shapefile and define extent of the figure

```{r, eval = FALSE}
# do intersection on data to limit kontur to portugal
st_pt <- st_intersection(data, pt)

# define aspect ratio based on bounding box
bb <- st_bbox(st_pt)

bottom_left <- st_point(c(bb[["xmin"]], bb[["ymin"]])) |> 
  st_sfc(crs = st_crs(data))

bottom_right <- st_point(c(bb[["xmax"]], bb[["ymin"]])) |> 
  st_sfc(crs = st_crs(data))

width <- st_distance(bottom_left, bottom_right)

top_left <- st_point(c(bb[["xmin"]], bb[["ymax"]])) |> 
  st_sfc(crs = st_crs(data))

height <- st_distance(bottom_left, top_left)

# handle conditions of width or height being the longer side
if (width > height) {
  w_ratio <- 1
  h_ratio <- height / width
} else {
  h_ration <- 1
  w_ratio <- width / height
}
```

### Convert it to a raster

```{r, eval = FALSE}
# convert to raster so we can then convert to matrix
size <- 5000

pt_rast <- st_rasterize(st_pt, 
                             nx = floor(size * w_ratio),
                             ny = floor(size * h_ratio))

mat <- matrix(pt_rast$population, 
              nrow = floor(size * w_ratio),
              ncol = floor(size * h_ratio))


```

### Create color palette and plot figure in high quality

```{r, eval = FALSE}
c1 <- met.brewer("OKeeffe2")
swatchplot(c1)

texture <- grDevices::colorRampPalette(c1, bias = 2)(256)
swatchplot(texture)

# plot the 3d map
rgl::close3d()

mat |> 
  height_shade(texture = texture) |> 
  plot_3d(heightmap = mat,
          zscale = 100/5,
          solid = FALSE, 
          shadowdepth = 0, #baseshape="hex", 
          theta = -10, phi = 60, zoom = .8)

outfile <- here('outputs/rayshader_PTplot.png')

{
  start_time <- Sys.time()
  cat(crayon::cyan(start_time), "\n")
  if (!file.exists(outfile)) {
    png::writePNG(matrix(1), target = outfile)
  }
  render_highquality(
    filename = outfile,
    interactive = FALSE,
    lightdirection = 280,
    lightaltitude = c(20, 80),
    lightcolor = c(c1[2], "white"),
    lightintensity = c(600, 100),
    samples = 450,
    width = 6000,
    height = 6000
  )
  end_time <- Sys.time()
  diff <- end_time - start_time
  cat(crayon::cyan(diff), "\n")
}
```

### Annotate the final figure

```{r, eval = FALSE}
img <- image_read(outfile)

colors <- met.brewer("OKeeffe2")
swatchplot(colors)

text_color <- darken(colors[7], .25)
swatchplot(text_color)

annot <- str_glue("This map shows population density of Portugal. ",
              "Population estimates are bucketed into 400 meter ",
              "squares.") |> 
  str_wrap(45)

img |> 
  image_crop(gravity = "center",
             geometry = "6000x3500+0-150") |> 
  image_annotate("Portugal Population Density",
                 gravity = "northwest",
                 location = "+200+100",
                 color = text_color,
                 size = 200,
                 weight = 700,
                 font = "Exo") |> 
  image_annotate(annot,
                  gravity = "west",
                  location = "+800+800",
                  color = text_color,
                  size = 125,
                  font = "Exo") |> 
  image_annotate(str_glue("Graphic by Ana Silva | ",
                      "Data: Kontur Population "),
                 gravity = "south",
                 location = "+0+100",
                 font = "Exo",
                 color = alpha(text_color, .5),
                 size = 70) |> 
  image_write(here('outputs/titled_rayshader_PTplot.png'))
```

![](outputs/titled_rayshader_PTplot.png)