Using ggmap and gganimation to visualize oil spill in Brazil coastline

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40

library(tidyverse)  # of course
library(rvest)      # to handle html scrapping
library(lubridate)  # handle datetime formats
library(glue)       # easily concat texts
library(knitr)      # markdown output tables
library(kableExtra) # formate markdown tables

# scrap and build url table
base_url <- "http://www.ibama.gov.br"

# page with update table
page <- read_html("http://www.ibama.gov.br/manchasdeoleo-localidades-atingidas")

# get the table
updates <- page %>% 
  html_table() %>% 
  .[[1]]

# get the href links
doc_links <- page %>%
  html_nodes(xpath = "//td/a") %>% 
  html_attr("href")

# put then togheter
doc_list <- updates %>% 
  set_names(c("date","description","items")) %>% 
  mutate( type = str_extract(items, "^[\\w]+") ) %>% 
  mutate( type = ifelse(type=="XLS","XLSX", type)) %>% # one file with "xls" extension
  mutate( link = paste0(base_url, doc_links) ) %>% 
  mutate( date=dmy(date) ) %>%
  as_tibble()

# save it, just in case we want to recover later
saveRDS(doc_list,"./data/oil_leakage_doc_list.rds")

# let's see the links list
doc_list %>% 
  head(10) %>% 
  kable() %>% 
  kableExtra::kable_styling(font_size = 10)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31

# check the pre existence of destination folders
XLSX_DEST_FOLDER <- "./data/oil_leakage_raw"
if(!file.exists("./data")) dir.create("./data")
if(!file.exists(XLSX_DEST_FOLDER)) dir.create(XLSX_DEST_FOLDER)

# download xlsx files
xlsx_filenames <- doc_list %>% 
  filter(type=="XLSX") %>% 
  select(date, link) %>% 
  arrange(date) %>% 
  split(1:nrow(.)) %>% 
  map_chr(function(.x){
    filename <- paste0(XLSX_DEST_FOLDER, "/",as.character(.x$date[1]), ".xlsx")
    download.file(url=.x$link[1], destfile = filename, mode = "wb")
    return(filename)
  })

# save an excel filename index
xlsx_index <- doc_list %>% 
  filter( type=="XLSX" ) %>% 
  arrange(date) %>% 
  mutate( filename = xlsx_filenames )

# save it
saveRDS(xlsx_index, "./data/excel_file_index.rds")

# let's see what we have
xlsx_index %>% 
  head(10) %>% 
  kable() %>% 
  kable_styling(font_size = 9)

1

library(xlsx) # import excel files (requires java/rjava)

1

## Warning: package 'xlsx' was built under R version 4.0.5

1
2
3
4
5
6

xlsx_file <- read.xlsx(file = "./data/2019-12-14.xlsx",sheetIndex = 1)

xlsx_file %>% 
  head(10) %>% 
  kable() %>% 
  kable_styling(font_size = 10)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12

# import each xlsx_file
files <- xlsx_index %>%
  split(1:nrow(.)) %>% 
  map(function(.x){
    print(.x$filename)
    read.xlsx(file = .x$filename[1], sheetIndex = 1, colClasses="character", stringsAsFactors=F) %>% 
      as_tibble() %>% 
      mutate(file.date=.x$date[1])
  })

# save it as cache
saveRDS(files, "./data/oil_leakage_imported_raw_excel.rds")

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135

# auxiliary function to translate the degree coordinates to decimal coordinates
translateCoord <- function(.coords){
  
  # extract orientation (N,S,W,E)
  directions <- str_extract(.coords,"\\w$")
  
  # split the numbers and rebuilds as xxDyy'zz"
  numbers <- .coords %>% 
    str_extract_all("(\\d\\.*)+") %>% 
    purrr::map(~paste0(.x[1],"d",.x[2],"'",.x[3],"\""))
  
  # modify char2dms to a safity call
  safe_char2dms <- safely(sp::char2dms, otherwise = NA)
  
  # readd orientation(NSWE) and converts to decimal
  numbers %>% 
    purrr::map2(directions, ~paste0(.x, .y)) %>% 
    purrr::map(safe_char2dms) %>% 
    purrr::map(purrr::pluck("result")) %>% 
    purrr::map(as.numeric) %>% 
    unlist() %>% 
    return()
}

# clean ant transform status data in to ordenated factor
statusToFactor <- function(.status){
  tibble( from_status = tolower(.status) )  %>% 
    mutate(
      to_status = case_when(
        str_detect(from_status, "manchas")   ~ "stains",
        str_detect(from_status, "esparsos")  ~ "traces/sparse",
        str_detect(from_status, "observado") ~ "not observed",
        T ~ as.character(NA)
      )
    ) %>% 
    mutate( fct_status = factor(to_status, levels = c("not observed", "traces/sparse", "stains"), ordered = T) ) %>% 
    pull(fct_status) %>% 
    return()
}

# read the raw excel importation
files <- readRDS("./data/oil_leakage_imported_raw_excel.rds")

# put all the excels in one dataframe
# perform a data cleanup
stains_raw <- files %>% 
  # the file form day 2019-11-06 hasn't coordinates
  keep(function(.x){ .x$file.date[1] != "2019-11-06"} ) %>% 
  # for each "excel" data
  map_df(function(.x){
    
    # some excels have two "municipio" columns
    if("municipio_" %in% names(.x)) .x <- select(.x, -municipio)
    
    # clean and normalize colnames
    col.names <- .x %>%
      names() %>% 
      tolower() %>%
      str_replace_all("munic.+pi.*", "municipio") %>% # acento no municipio 
      str_replace_all("name", "localidade") %>% 
      str_replace_all("latitutde","latitude")
    
    # rename the cols and select what we'll use
    .x %>% 
      set_names(col.names) %>% 
      select(file.date, localidade, municipio, 
             estado, data_avist, data_revis,
             latitude, longitude, status) %>% 
      
      # handle portuguese accents
      mutate(
        localidade = iconv(localidade, from="UTF-8", to="LATIN1"),
        municipio  = iconv(municipio, from="UTF-8", to="LATIN1"),
        estado     = iconv(estado, from="UTF-8", to="LATIN1"),
        status     = iconv(status, from="UTF-8", to="LATIN1"),
        data_avist = ymd(as.character(data_avist)), 
        data_revis = ymd(as.character(data_revis))
      ) %>% 
      return()
  })

# we have to translate degree coordinates to decimal
stains_coord <- stains_raw %>%   
  mutate(
    lat.degree = translateCoord(latitude),
    lon.degree = translateCoord(longitude)
  ) %>% 
  # in some of the files, the coordinates are indeed in decimal 
  mutate(
    lat.degree = ifelse(is.na(lat.degree), as.numeric(latitude), lat.degree),
    lon.degree = ifelse(is.na(lon.degree), as.numeric(longitude), lon.degree)
  )

# in some files, the "UF" column has the name of the states 
# and others the code of states, we handle this making a name<->code table
estado_uf <- files %>% 
  tail(1) %>% 
  .[[1]] %>% 
  mutate(
    estado = iconv(estado, from="UTF-8", to="LATIN1")
  ) %>% 
  select(estado, sigla_uf) %>% 
  distinct() %>% 
  arrange(estado)

# join the state name<->code table
stains_clean <- stains_coord %>% 
  left_join(estado_uf, by = "estado") %>% 
  # handle when the info is correct/incorrect
  mutate(sigla_uf = ifelse(is.na(sigla_uf),estado,sigla_uf)) %>% 
  # there is files without state code
  filter(!is.na(sigla_uf)) %>% 
  select(-estado) %>% 
  # rebuild state name
  inner_join(estado_uf, by="sigla_uf") %>% 
  # reselect only interesting cols
  select(file.date, data_avist, data_revis, sigla_uf, estado, municipio, localidade, lat.degree, lon.degree, status)

# transform the "status" info in a factor column
stains <- stains_clean %>% 
  mutate( status = statusToFactor(status) ) %>% 
  filter(complete.cases(.)) %>% 
  # english col names
  set_names(c("file.date", "sighting_date","revision_date","uf","state","county",
              "local","lat.degree","lon.degree","status")) %>% 
  # last clean-up
  filter(
    lon.degree < 0, # anything greater than that is a error
    revision_date < ymd(20200101), # anything greater than that is a error
    sighting_date < ymd(20200101), # anything greater than that is a error
    complete.cases(.) # missing data is not relevant
  )

# save it, just in case
saveRDS(stains, "./data/oil_leakage.rds")

1
2
3
4
5
6
7

stains <- readRDS("./data/oil_leakage.rds")

# checking the final format
stains %>% 
  head(10) %>% 
  kable() %>% 
  kable_styling(font_size = 8)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21

# ggmap
library(ggmap)
library(yaml) # used to not version the google's map key

# before use ggmap with google maps it's necessary
# read and register a Key for Google Map API
# config <- yaml::read_yaml("./config.yml")
# register_google(key=config$google_map_key)

# get the map area
bbox <- make_bbox(lon = stains$lon.degree, lat=stains$lat.degree, f = .1)
gmap <- get_map(location=bbox, source="google", maptype="terrain")

# plot the data of observations on revision date 2019-11-10
# of location with stains or sparse residuos
ggmap(gmap) +
  geom_point(data=filter(stains, status!="not observed", revision_date==ymd(20191110)),
             aes(x=lon.degree, y=lat.degree, color=status), size=2, alpha=.8) +
  scale_color_manual(values=c("#888888","black","green")) +
  theme_void() +
  ggplot2::coord_fixed() 

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20

library(gganimate)
bbox <- make_bbox(lon = stains$lon.degree, lat=stains$lat.degree, f = .1)
gmap <- get_map(location=bbox, source="google", maptype="terrain")

p <- ggmap(gmap) +
  geom_point(data=filter(stains, status!="not observed"),
             aes(x=lon.degree, y=lat.degree,
                 color=status, group=county),
             size=3, alpha=.8) +
  scale_color_manual(values=c("#888888","black")) +
  theme_void()

anim <- p  +
  theme( legend.position = "bottom" ) +
  transition_time(revision_date) +
  labs(title="Oil Sighting on Brazil Coast",
       subtitle = "Date:  {frame_time} - source: IBAMA") +
  shadow_wake(wake_length = .15)

anim

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19

# let's organize the "county/locality" by latitude
stains %>% 
  select(local, lat.degree, revision_date, status) %>% 
  group_by(local) %>% 
  # calculates a unique latitude for local
  mutate( avg.lat = mean(lat.degree) ) %>% 
  ungroup() %>% 
  # put the locals in latitude order
  mutate( local = fct_reorder(local, avg.lat) ) %>% 
  select( local, revision_date, status ) %>% 
  distinct() %>% 
  # plot as a tile chart (local x date)
  arrange( local, revision_date ) %>% 
  ggplot(aes(x=revision_date, y=local)) +
  geom_tile(aes(fill=status)) +
  scale_fill_manual(values = c("lightgreen", "#888888", "black")) +
  theme_minimal() +
  theme( axis.text.y = element_blank() ) +
  labs(title="Oil Sighting by Latitude and Date", x="date", y="local/latitude")