Animations

libraries:

library(tidyverse)
library(ggthemes) 
library(plotly)
library(sf)
library(leaflet)
library(gganimate)
library(tmap)

Animation is a method in which pictures are manipulated to appear as moving images (link). Basically it’s a sequence of images (rapid succession of sequential images).

Actually we already tried different animations during the previous practical sessions. But let’s look it now more closely. In general animations can be created for two purposes:

• process animation (changes in time),
• animation of static plot to reveal different angles.

Animation as sequence of images

In next example we try to create the animation, where we demonstrate how the contour of Estonia is changing during the simplification process:

Download the Estonian contour:

download.file("http://aasa.ut.ee/GisMaps/data/est_wgs84.zip", destfile = "est_wgs84.zip") 
unzip("est_wgs84.zip")

est_cont <- st_read("eestimaa_wgs84.shp")

## Reading layer `eestimaa_wgs84' from data source 
##   `C:\ANTO\loengud\geopythonR\rspatial_Git\rspatial\eestimaa_wgs84.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 1 field
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 21.76431 ymin: 57.50931 xmax: 28.209 ymax: 59.82202
## Geodetic CRS:  WGS 84

# transform it to 3301:
est_cont_3301 <- st_transform(est_cont, 3301)

# simplfy /generalize:
est_cont_3301_smpl <- st_simplify(est_cont_3301, preserveTopology = F, dTolerance = 200)

Now your task is to create sequence of images, where every next image is a bit more simplified than previous. First we create the vector of simplification level (dTolerance) for every step :

# create data frame, where first value of dTolerance is 102
steps <- data.frame(step = 102) # 102 comes after testinfg several different values

# Run the for/next loop to create the vector:
for(j in 1 : 100){
  steps_tmp <- steps[j,] # take last dTolerance value
  steps_tmp <- 1.07 * steps_tmp # multiply the value by 1.02 (+2%)
  steps <- rbind(steps, data.frame(step = steps_tmp)) # bind the new value to data frame as new row
}

First 10 steps:

head(steps, n = 10)

##        step
## 1  102.0000
## 2  109.1400
## 3  116.7798
## 4  124.9544
## 5  133.7012
## 6  143.0603
## 7  153.0745
## 8  163.7897
## 9  175.2550
## 10 187.5228

Last 10 steps:

tail(steps, n = 10)

##         step
## 92  48141.98
## 93  51511.92
## 94  55117.75
## 95  58975.99
## 96  63104.31
## 97  67521.62
## 98  72248.13
## 99  77305.50
## 100 82716.88
## 101 88507.07

In the next step we create a fixed frame around Estonia. Otherwise, the animation will start “jumping around”. For that we use a bounding box which will be converted to a polygon:

a <- st_bbox(est_cont_3301_smpl)

# Result: 
a

##      xmin      ymin      xmax      ymax 
##  369047.8 6377146.4  739155.1 6634019.1

# convert it to dataframe:
a <- data.frame(x = c(a[1], a[3]), y = c(a[2], a[4]))
a

##             x       y
## xmin 369047.8 6377146
## xmax 739155.1 6634019

# sf object:
a <- st_as_sf(a, coords=c("x", "y"), crs = 3301)

Plot the bbox and contour of Estonia:

tmap_mode("plot") # confirm the static mode

## tmap mode set to plotting

tm_shape(a)+
  tm_dots("red", size= 2, shape = 3)+
  tm_shape(est_cont)+
  tm_borders()

To keep stuff in folder tidy, we create sub-directory for animation frames:

dir.create("./anim")

In next stage we produce animation frames. We are using the simple for/next loop, where simplified contour is created for every simplification step:

# create the image/map for every generalization step:
for(i in 1 : nrow(steps)){
  est_cont_ee <- est_cont_3301 # create copy of contour
  est_cont_ee <- st_simplify(est_cont_ee, dTolerance = steps[i,1]) # simplification for level i
  
  # name of image (number in sequence, based on this numbering the images are bound together):
  frameName <- str_pad(i, 4, side = "left", pad = "0")
  
  
  # previous step is needed because numbering in file name should be in format 0008, ,009, 0010 etc (not: 8, 9, 10 etc)!
    myplot <- ggplot()+
      theme_minimal()+
      geom_sf(data = a, colour="white")+
      geom_sf(data = est_cont_ee)
  
  # save the frame
  ggsave(filename = paste0("./anim/anim_", frameName, ".png"), plot = myplot, dpi = 100, width = 5, height = 4, units = "in")  
  
  # visualize the progress:
  print(i)
  
}

Now you should have 101 images of Estonian contour.

Next stage assumes, that you have ffmpeg installed in your computer and it’s accessible from command line. Step-by-step tutorial How to install FFmpeg on Windows.
ffmpeg doesn’t have the graphical user interface (GUI). It can be used from terminal (e.g. Command Prompt), but you can run it directly from R as well:

#Terminal:
cd C:/ANTO/loengud/geopythonR/2019/rmd/anim/
ffmpeg -start_number 1 -i anim_%04d.png -vcodec libx264 -crf 25  -pix_fmt yuv420p test.mp4


#RStudio: 
system("cd C:/ANTO/loengud/geopythonR/2019/rmd/anim/")
system("ffmpeg -start_number 1 -i anim_%04d.png -vcodec libx264 -crf 25  -pix_fmt yuv420p test.mp4")

FFmpeg is very powerful software. You can create animations, add soundtrack, titles etc. Unfortunately it’s not intuitive and first steps can be really hard. But if you mange, you can create something like that:

Result in Vimeo: link

If ffmpeg is for the beginner a bit hard, then magick can help as well. Here is an example how to produce GIF:

library(magick)

## Warning: package 'magick' was built under R version 4.1.1

## Linking to ImageMagick 6.9.12.3
## Enabled features: cairo, freetype, fftw, ghostscript, heic, lcms, pango, raw, rsvg, webp
## Disabled features: fontconfig, x11

## list file names and read in
imgs <- list.files("./anim/", full.names = TRUE)
img_list <- lapply(imgs, image_read)

## join the images together
img_joined <- image_join(img_list)

## animate at 2 frames per second
img_animated <- image_animate(img_joined, fps = 10)

## view animated image
img_animated

## save to disk
image_write(image = img_animated,
            path = "anim.gif",)

Let’s leave it for now and look some other examples.

Sinking Estonia

In next example we try to animate the change of Estonian contour caused by rise of sea level. Example itself is quite primitive again and unrealistic, but good enough to learn.
For the animation we need the elevation layer of Estonia. With help of packages called raster we can easily download it directly from web:

library(raster)
est_elev <- getData('alt', country = 'EST', mask = F) 
plot(est_elev)

class(est_elev)

## [1] "RasterLayer"
## attr(,"package")
## [1] "raster"

In next step we create animation steps. It means we create the sequence of height above sea level in meters (from 0 to 280 meters)

water_step <- seq(0, 280, 10) # after every 10 meters

water_step

##  [1]   0  10  20  30  40  50  60  70  80  90 100 110 120 130 140 150 160 170 180
## [20] 190 200 210 220 230 240 250 260 270 280

Every value in this vector corresponds to the possible sea level rise in meters.
Now design the map (again just for one frame is enough):

# calculate the elevation after sea level rise ( water_step[10] corresponds to 90 meters):
  est_elev_tmp <- est_elev - water_step[10]
  
  # create masking layer:
  est_elev_mask <- est_elev %in% (water_step[10]:300) 

Quick glimpse:

plot(est_elev_mask)

Thats the elavation map where the elavation is more than 90 m.a.s.l. [meters above the sea level]
Design the single animation frame (> 90 masl):

  tm_shape(est_elev_tmp)+
    tm_raster(style  = "fixed", 
            breaks = c(20, 50, 75, 100, 150, 200, 250), 
            palette = tmaptools::get_brewer_pal("-RdYlGn", n = 7, contrast = c(0, 1)),
            title = "Estonian elevation, (m.a.s.l.)", 
            midpoint = NA)+
  tm_shape(est_elev_mask)+
    tm_raster(alpha = c(1, 0), legend.show = F)+
  tm_shape(est_cont)+
    tm_borders()+
  tm_credits(paste0("Estonian countour after sea level rise of ", water_step[10] ," meters"), position=c("left", "bottom"))

Just to be sure, that tmap is in static mode:

tmap_mode("plot")

## tmap mode set to plotting

for(i in 1:length(water_step)){
  # calculate the elevation after sea level rise: 
  est_elev_tmp <- est_elev - water_step[i]
  
  # create masking layer:
  est_elev_mask <- est_elev %in% (water_step[i]:300)
  plot(est_elev_mask)
  
  # numbering for animation frame image:
  j <- ifelse(nchar(i) == 1, paste0("0", i), as.character(i))
  file_name <- paste0("est_elev_", j, ".png")
  
  # create the single animation frame:
  map <- tm_shape(est_elev_tmp)+
    tm_raster(style  = "fixed", 
            breaks = c(20, 50, 75, 100, 150, 200, 250), 
            palette = tmaptools::get_brewer_pal("-RdYlGn", n = 7, contrast = c(0, 1)),
            title = "Estonian elevation, (m.a.s.l.)", 
            midpoint = NA)+
  tm_shape(est_elev_mask)+
    tm_raster(alpha = c(1, 0), legend.show = F)+
  tm_shape(est_cont)+
    tm_borders()+
  tm_credits(paste0("Estonian countour after sea level rise of ", water_step[i] ," meters"), position=c("left", "bottom"))
  
  # save it:
  tmap_save(map, filename = file_name)

}

In next stage we indicate to animation frames and bind them together into one GIF-animation:

library(purrr)
library(magick)

list.files(pattern = "est_elev_") %>%  # list of animation frames
  map(image_read) %>% # reads each path file
  image_join() %>% # joins image
  image_animate(fps=5) %>% # animates
  image_write("est_sinking.gif") # write to current directory

And the results should be something like that:

Animate movement

The example is very abstract again. Let’s assume that we have the data about unkown flying object which is making circles over Estonia. In fact you create the dataset. First calculate the centroid of Estonia and make 100 km buffer around it:

est_cntr_sf <- st_centroid(est_cont_3301_smpl)
st_crs(est_cntr_sf)

## Coordinate Reference System:
##   User input: EPSG:3301 
##   wkt:
## PROJCRS["Estonian Coordinate System of 1997",
##     BASEGEOGCRS["EST97",
##         DATUM["Estonia 1997",
##             ELLIPSOID["GRS 1980",6378137,298.257222101,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4180]],
##     CONVERSION["Estonian National Grid",
##         METHOD["Lambert Conic Conformal (2SP)",
##             ID["EPSG",9802]],
##         PARAMETER["Latitude of false origin",57.5175539305556,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8821]],
##         PARAMETER["Longitude of false origin",24,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8822]],
##         PARAMETER["Latitude of 1st standard parallel",59.3333333333333,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8823]],
##         PARAMETER["Latitude of 2nd standard parallel",58,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8824]],
##         PARAMETER["Easting at false origin",500000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8826]],
##         PARAMETER["Northing at false origin",6375000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8827]]],
##     CS[Cartesian,2],
##         AXIS["northing (X)",north,
##             ORDER[1],
##             LENGTHUNIT["metre",1]],
##         AXIS["easting (Y)",east,
##             ORDER[2],
##             LENGTHUNIT["metre",1]],
##     USAGE[
##         SCOPE["Topographic mapping (large scale)."],
##         AREA["Estonia - onshore and offshore."],
##         BBOX[57.52,20.37,60,28.2]],
##     ID["EPSG",3301]]

est_cntr_buf <- st_buffer(est_cntr_sf, dist = 100000)

Plot it:

tm_shape(est_cont_3301_smpl)+
  tm_borders()+
tm_shape(est_cntr_buf)+
  tm_borders("red", lwd = 2)

The circle is a polygon defined by points which are connected with each other. Convert the polygon to points:

# get coordinates as data frame:
est_cntr_buf_df <- st_coordinates(est_cntr_buf) %>% 
  as_tibble()

# rank the coordinates:
est_cntr_buf_df <- est_cntr_buf_df %>% 
  mutate(rnk = seq(1, nrow(est_cntr_buf_df), 1))

# plot it:
ggplot()+
  theme_minimal()+
  geom_sf(data= est_cont_3301_smpl, col = "grey", fill = "grey97", size = 0.25)+
  geom_point(data = est_cntr_buf_df, aes(x = X, y = Y, size= rnk, alpha= rnk), col = "red")

As you notice we are using alpha to add the fading tail to our moving dots. Currently it’s too long. We should add it to 10 location points. Now our task is to create animation frames. But firstly we should design the single animation frame:

# first 10 points:
est_cntr_buf_df_tmp <- est_cntr_buf_df[1:10,]

# firste frame:
ggplot()+
  theme_minimal()+
  geom_sf(data= est_cont_3301_smpl, col = "grey", fill = "grey97", size = 0.25)+
  geom_path(data = est_cntr_buf_df, aes(x = X, y = Y), col = "pink", size= 0.1)+
  geom_point(data = est_cntr_buf_df_tmp, aes(x= X, y = Y, alpha = rnk, size= rnk), colour = "red")+
  scale_alpha(range = c(0, 1))+
  scale_size(range = c(0.5, 4))

Our dataset contains 121 location points. This means we can create single animation frame for every dot. But firstly we have to create again the numbering of images:

### calculate the number of characters in rank (new column):
est_cntr_buf_df <- est_cntr_buf_df %>% 
  mutate(rnk_nchar = nchar(rnk))

# create separate data frame for image name prefixes:
rnk_joiner <- data.frame(rnk_nchar = c(1, 2, 3),
                         prefix = c("00", "0", NA))

rnk_joiner

##   rnk_nchar prefix
## 1         1     00
## 2         2      0
## 3         3   <NA>

## join two table:
est_cntr_buf_df <- left_join(est_cntr_buf_df, rnk_joiner)

## Joining, by = "rnk_nchar"

head(est_cntr_buf_df)

## # A tibble: 6 x 7
##         X        Y    L1    L2   rnk rnk_nchar prefix
##     <dbl>    <dbl> <dbl> <dbl> <dbl>     <dbl> <chr> 
## 1 685719. 6504580.     1     1     1         1 00    
## 2 685582. 6499346.     1     1     2         1 00    
## 3 685171. 6494127.     1     1     3         1 00    
## 4 684488. 6488936.     1     1     4         1 00    
## 5 683534. 6483789.     1     1     5         1 00    
## 6 682312. 6478698.     1     1     6         1 00

# calculate the rank as character:
est_cntr_buf_df <- est_cntr_buf_df %>% 
  mutate(frame_name = ifelse(!is.na(prefix), paste0(prefix, rnk), rnk))

Create animation frames and save them as png-files:

for(i in 1: nrow(est_cntr_buf_df)){
  est_cntr_buf_df_tmp2 <- est_cntr_buf_df %>% 
    filter(rnk > i - 1 & rnk < i +10) 
  
  est_cntr_buf_df_head <- est_cntr_buf_df %>% 
    filter(rnk == i + 10)
  
  ggplot()+
    theme_light()+
    geom_sf(data= est_cont_3301_smpl, col = "grey", fill = "grey97", size = 0.25)+
    geom_path(data = est_cntr_buf_df, aes(x = X, y = Y), col = "pink", size= 0.1)+
    geom_point(data = est_cntr_buf_df_tmp2, aes(x= X, y = Y, alpha = rnk, size= rnk), colour = "orange")+
    geom_point(data = est_cntr_buf_df_head, aes(x= X, y = Y), colour = "darkred", size= 2.5)+
    scale_alpha(range = c(0, 0.7))+
    scale_size(range = c(0.2, 2))+
    guides(alpha = F, size = F)
  
  frame_name <- paste0("airplane_", est_cntr_buf_df$frame_name[i], ".png")
  ggsave(filename = frame_name, dpi = 200, height = 3, width = 4.5) 
}

For animatin we use some magick:

Bind the images together as animation:

getwd()

list.files(pattern = "airplane_") %>% 
  map(image_read) %>% # reads each path file
  image_join() %>% # joins image
  image_animate(fps = 25) %>% # animates ('fps' defeines the speed of animation, frames per second)
  image_write("est_UFO.gif") # write to current dir

Same method is applied in here:

ggplot2 & rayshader

Elevated population map

During the previous session the rayshader library was introduced. It’s very powerful and showy (the installation inferno is worth it).

In next example we look, how to convert “regular” plane ggplot2 plot to 3D.

First, download the data if needed:

download.file("http://aasa.ut.ee/GisMaps/data/population_2017.zip", destfile = "population_2017.zip") 
unzip("population_2017.zip")

Import it:

popul_2017 <-  st_read("population_2017.shp", quiet = T)
st_crs(popul_2017) <- 3301

# calculate population ranks for the municipalities:
popul_2017 <- popul_2017 %>% 
  mutate(pop_rank = min_rank(desc(popultn))) 

# First you have to convert the geolayer from Estonian CRS (**LEST97, epsg:3301**) to global geographical (**WGS84; epsg:4326**) coordinate system:
popul_2017_wgs84 <- st_transform(popul_2017, 4326)


# check the structure:
glimpse(popul_2017_wgs84)

## Rows: 213
## Columns: 9
## $ OKOOD    <chr> "0728", "0490", "0580", "0868", "0303", "0689", "0478", "0349~
## $ county   <chr> "Harju", "Viljandi", "Harju", "Harju", "Pärnu", "Saare", "Saa~
## $ mncplty  <chr> "Saue", "Mõisaküla", "Paldiski", "Vasalemma", "Kihnu", "Ruhnu~
## $ popultn  <dbl> 5809, 784, 3806, 2499, 701, 147, 1887, 13635, 1817, 876, 706,~
## $ voters   <dbl> 4497, 715, 3055, 1994, 624, 142, 1687, 11325, 1625, 748, 620,~
## $ munc_cd  <dbl> 728, 490, 580, 868, 303, 689, 478, 349, 550, 634, 386, 721, 8~
## $ cont_cd  <dbl> 37, 84, 37, 37, 67, 74, 74, 74, 74, 74, 74, 74, 74, 37, 37, 5~
## $ geometry <MULTIPOLYGON [arc_degree]> MULTIPOLYGON (((24.54087 59..., MULTIPO~
## $ pop_rank <int> 36, 188, 62, 82, 197, 212, 107, 13, 109, 179, 196, 155, 211, ~

Draw the map:

map1 <- ggplot(data = popul_2017_wgs84) +
  geom_sf(aes(fill = popultn), color = NA, size = 0.3) +
  scale_fill_viridis_c(option = "inferno",
                    name = "Population",
                    
                    trans= 'log10') + 
  theme_void()+
  guides(fill = guide_colourbar(title.position = "top"))+
  ggtitle("Population size in 2017") 

map1

Avoid scientific notation:

options(scipen = 999)

map1

And now we can convert the map to 3D model, where municipalities are elevated according to the population size:

rayshader::plot_gg(map1,width = 3.5,windowsize = c(1000, 800), multicore = TRUE,
        zoom = 0.8, phi = 35, theta = 30, sunangle = 225, soliddepth = -100)

rayshader::render_snapshot()

Render it to the rotating animation:

rayshader::render_movie("movie_tmp.mp4",frames = 720, fps=30,zoom=0.6,fov = 30)

It should be in your working directory (getwd()) now.

3D elevation model of Estonia

Start additional libraries:

library(RColorBrewer)
library(raster)

Get elevation model and plot it:

est_elev <- getData('alt', country = 'EST', mask = F) 

plot(est_elev)

Convert it to data frame (XYZ-data):

# Convert raster to a DataFrame:
est_elev_df <- rasterToPoints(est_elev, spatial = F) %>% 
  as_tibble()

glimpse(est_elev_df)

## Rows: 143,652
## Columns: 3
## $ x       <dbl> 28.12083, 28.10417, 28.11250, 28.12083, 28.12917, 28.13750, 28~
## $ y       <dbl> 59.79583, 59.78750, 59.78750, 59.78750, 59.78750, 59.78750, 59~
## $ EST_alt <dbl> 0, 6, 6, 2, 3, 3, 2, 0, 1, 0, 9, 17, 9, 6, 10, 20, 16, 17, 17,~

Create the colour palette with 11 color steps:

myPalette <- colorRampPalette(rev(brewer.pal(11, "RdYlGn")))  # Create color palette

Plot the map:

map2 <- ggplot()+
  geom_raster(data = est_elev_df,aes(x = x, y = y, fill = EST_alt))+
  #geom_sf(data=est_cont, fill=NA, color="black", size = 0.25)+
  scale_fill_gradientn(colours = myPalette(4), name = "Elevation (m)")+
  theme_void()+
  theme(plot.margin = grid::unit(c(0,0,0,0), "mm"),
        legend.position = "bottom",
        legend.text = element_text(family = "mono"),
        legend.title = element_text(family = "mono"))

map2

Convert it to 3D:

rayshader::plot_gg(map2,
        #width = 3.5,
        windowsize = c(1000, 1000), 
        shadow_intensity=0.4,
        multicore = TRUE,
        zoom = 0.6)

rayshader::render_snapshot()

Render it to the rotating animation:

rayshader::render_movie("movie_estElev.mp4",frames = 720, fps = 30, zoom = 0.6,fov = 30)

It should be in your working directory (getwd()) now.

---

Author: Anto Aasa, Iuliia Burdun
Supervisors: Anto Aasa & Lika Zhvania
LTOM.02.041
Last update: 2022-12-07 11:09:52

.