How to make the 3-D plots shown here: https://twitter.com/grant_mcdermott/status/1149425928176431104

gfw2016-rayshader.R

### LOAD PACKAGES ###

# install.packages("pacman")
pacman::p_load(tidyverse, bigrquery, DBI, sf, stars, rayshader)


### DOWNLOAD GLOBAL FISHING WATCH DATA FROM GOOGLE BIGQUERY ###

## You'll need a valid BigQuery billing ID for this next part. Luckily, you can
## easily sign up for a 12-month free trial that gives you access to BQ, as well
## as all of the other products on Google Cloud Platform. Moreover, you get 1 TB
## of free queries every month, regardless of whether your free trial has expired
## or not. I'm going to skip over most other details of this section. But for
## details on how to sign up for the free 12-month trial, as well as more 
## information about what I'm doing and why, see my databases lecture (#16):
## https://github.com/uo-ec607/lectures#lecture-outline-and-quicklinks

# library(tidyverse) ## Already loaded
# library(bigrquery) ## Already loaded
# library(DBI) ## Already loaded

## Connect to GFW data on BigQuery

billing_id = Sys.getenv("GCE_DEFAULT_PROJECT_ID") ## Replace with your project ID

gfw_con = 
  dbConnect(
    bigrquery::bigquery(),
    project = "global-fishing-watch",
    dataset = "global_footprint_of_fisheries",
    billing = billing_id
    )

## Reference the fishing_effort table
effort = tbl(gfw_con, "fishing_effort")

## Define the desired bin resolution in degrees
resolution = 1

## Download binned (i.e. gridded) data, summed over all of 2016 
gfw16 =
  effort %>%
  filter(
    `_PARTITIONTIME` >= "2016-01-01 00:00:00",
    `_PARTITIONTIME` <= "2016-12-31 00:00:00"
    ) %>%
  filter(fishing_hours > 0) %>%
  mutate(
    lat_bin = lat_bin/100,
    lon_bin = lon_bin/100
    ) %>%
  mutate(
    lat_bin_center = floor(lat_bin/resolution)*resolution + 0.5*resolution,
    lon_bin_center = floor(lon_bin/resolution)*resolution + 0.5*resolution
    ) %>%
  group_by(lat_bin_center, lon_bin_center) %>%
  summarise(fishing_hours = sum(fishing_hours, na.rm=T)) %>%
  collect()


### CONVERT TO SF (VECTORISED) OBJECT FOR PLOTTING ###

## To get the above (gridded) data frame into a nice sf object, we have to first
## go through the rigmoral of converting it into a raster object (using
## rasterFromXYZ), then a stars object... and then an sf object. I'll also go 
## ahead and set the CRS to the equal earth projection.

gfw16_sf = 
  gfw16 %>%
  raster::rasterFromXYZ(crs= "+proj=longlat +datum=WGS84 +no_defs") %>%
  st_as_stars() %>%
  st_as_sf() %>%
  st_transform("+proj=eqearth +wktext") 
  

## Generate and save a normal ggplot object
p = 
  gfw16_sf %>%
  mutate(log_fishing_hours = log(fishing_hours)) %>%
  ggplot() +
  geom_sf(aes(col=log_fishing_hours, fill=log_fishing_hours)) +
  scale_colour_viridis_c(name = "Log fishing hours")  +
  scale_fill_viridis_c(name = "Log fishing hours") +
  labs(
    title = "Global fishing effort in 2016",
    subtitle = paste0("Effort binned at the ", resolution, "° level."),
    y = NULL, x = NULL,
    caption = "Data from Global Fishing Watch"
  ) +
  theme_void()

## 2-D plot
p 

######################################
### ADD 3-D EFFECTS WITH RAYSHADER ###

# library(rayshader) ## Already loaded

## Generate interactive plot (keep this open)
plot_gg(p, zoom = 0.8, multicore = TRUE, width = 8 ,height = 4)

## Animated MP4 ##

## The `rayshader::render_movie()` defaults work pretty well, but I want to rotate, 
## tilt, and zoom in particular ways. So I'll specify these parameters separately.
num_frames = 360
phivec = 30 + 60 * 1/(1 + exp(seq(-7, 20, length.out=num_frames/2)/2))
phivec = c(phivec, rev(phivec))
thetavec = 60 * sin(seq(0, 359, length.out=num_frames) * pi/180)
zoomvec = 0.45 + 0.2 * 1/(1 + exp(seq(-5, 20, length.out=num_frames/2)))
zoomvec = c(zoomvec, rev(zoomvec))
## Render animation
render_movie(
  filename = "gfw2016.mp4", type = "custom", 
  phi = phivec, zoom = zoomvec, theta = thetavec,
  frames = num_frames,  fps = 30
  )
  
## Alternatively (or in addition), you could generate the frames for the movie
## separately as PNGs. Then you can convert to MP4 (or GIF or whatever) using
## ffmeg.
## Not run
# dir.create("frames") ## Make sure directory is created
# lapply(1:num_frames, function(i) {
#   render_camera(theta = thetavec[i],phi = phivec[i],zoom = zoomvec[i])
#   render_snapshot(paste0("./frames/frame", i, ".png"))
#   })
# ## Shell (bash) call to ffmeg
# system("ffmpeg -framerate 30 -i ./frames/frame%d.png -vcodec libx264 gfw2016.mp4")
## End not run


## Second plot: Not log transformed ##

pn = 
  gfw16_sf %>%
  # mutate(log_fishing_hours = log(fishing_hours)) %>%
  ggplot() +
  geom_sf(aes(col=fishing_hours, fill=fishing_hours)) +
  scale_colour_viridis_c(name = "Fishing hours")  +
  scale_fill_viridis_c(name = "Fishing hours") +
  labs(
    title = "Global fishing effort in 2016",
    subtitle = paste0("Effort binned at the ", resolution, "° level."),
    y = NULL, x = NULL,
    caption = "Data from Global Fishing Watch"
  ) +
  theme_void()

plot_gg(pn, zoom = 0.5, multicore = TRUE, width = 8 ,height = 4, fov = 70)

## Parameters
num_frames = 360
phivec = 90 * 1/(1 + exp(seq(-7, 20, length.out=num_frames/2)/2))
phivec = c(phivec, rev(phivec))
thetavec = 360 - 90 * sin(seq(0, 180, length.out=num_frames/3) * pi/180)
thetavec = c(rep(0, num_frames/3), thetavec, rep(0, num_frames/3))
zoomvec = 0.25 + 0.25 * 1/(1 + exp(seq(-5, 20, length.out=num_frames/2)))
zoomvec = c(zoomvec, rev(zoomvec))
## Render animation
render_movie(
  filename = "gfw2016-nolog.mp4", type = "custom", 
  phi = phivec, theta = thetavec, zoom = zoomvec,
  frames = num_frames,  fps = 30
)