Example code to generate animation frames of Himawari-8 hotspots

him_anim.r

# General data manipulation
library(tidyverse)

# Spatial data
library(sf)

# Thematic maps
library(tmap)

# Parallel processing
library(furrr)

# Date manipulation
library(lubridate)

# GIS Raster
library(raster)

# Enable parallel processing of output
plan(multiprocess)


############
## Hotspots are available from JAXA's Himawari-8
## You can request a research login and download
## Monthly and Daily data files from their FTP server
##
## https://www.eorc.jaxa.jp/ptree/registration_top.html
##
## Once you are registered, files can be downloaded from:
##
## ftp://ftp.ptree.jaxa.jp/pub/himawari/L3/WLF/bet
##
## Monthly files are only available at the end of each month
## For months in progress download the daily files instead

# After you have downloaded your files
# Load all monthly files from the monthly directory and merge  
fls = list.files("monthly/",full.names = TRUE)
all = map(fls,read_csv,quote="'")
data = bind_rows(all)

# Load all daily files from the daily directory and merge
fls = list.files("daily/",full.names = TRUE)
all = map(fls,read_csv,quote="'")
data2 = bind_rows(all)

# Combine monthly and daily files - the format is the same
data = bind_rows(data,data2)

# Himawari-8 data is the entre disk at longitude 140.7E
# Here we define an extent for Australia to crop it to
ex = extent(c(xmin=-112,xmax=155,ymin=-44,ymax=-10))

# Convert data to a spatial object
data = st_as_sf(data,crs=4326,wkt="pixel")

# Crop to Australia
data = st_crop(data,ex)

# Simplify to fields of interest
data=dplyr::select(data,"#obstime",lon,lat,viewzenang,viewazang,pixwid,pixlen,fire_idx,sunglint,firepower,freq,fQC)
names(data)[1]="ObsTime"

# Generate fields for year, month, day, hour
data = data %>% mutate(year = year(ObsTime),day = day(ObsTime),month = month(ObsTime), hour=hour(ObsTime))

# Strip the geometry at this point, for simpler processing.
st_geometry(data) = NULL

# Create a time field to group the hotspots by hour, and then group by hour, lat and long
# And generate sumamry statistics (mean and maximum radiative power) for each coordinate and hour
data = data %>% 
  mutate(TS = make_datetime(year,month,day,hour,0)) %>% 
  group_by(TS,lat,lon) %>% summarise(cnt = length(firepower),mean = mean(firepower),max=max(firepower))

# Also generate and overall "all-time" summary while we're here, removing the time element
datax = data %>% 
  group_by(lat,lon) %>% 
  summarise(count = length(cnt),mean=mean(mean),max=max(max))

datax = st_as_sf(datax,coords=c("lon","lat"),crs=4326)



# Define the start and end dates we want to plot - this should a subset
# of the dates of the data you downloaded
st_date = as.POSIXct("2019-10-15")
en_date = max(data$TS)
data2 = filter(data,TS > st_date)
unq_times = seq(st_date,en_date,by=60*60)

# During some timesteps we have no hotspots - when this happens
# This is a "dummy" set of hotspots to plot, just to avoid messing up
# the plot algorithm 
ddata = tibble(x=0,y=0,firepower=0,max=0)
ddata = st_as_sf(ddata,coords=c("x","y"),crs=4326)

# This is a GeoTIFF background satellie image of the area of interest
syd = stack("tst.tiff")

# It happens to be a four-band TIFF with transparency - just keep the first
# three bands and discard transparency
syd = subset(syd,c(1,2,3))


# Main plotting function
proc=function(tt){
  # tdata = hotspots we will plot for this timestep
  tdata=filter(data2,TS == unq_times[tt])
  
  # fdata = all previous hotspots, which we will plot in grey 
  fdata = filter(data2, TS <= unq_times[tt])
  
  # If either of these are empty, use the blank dummy dataset instead
  if(nrow(tdata)==0){
    tdata = ddata
  }
  if(nrow(fdata)==0){
    fdata = ddata
  }
  
  # Make them spatial objects
  tdata = st_as_sf(tdata,coords=c("lon","lat"),crs=4326)
  fdata = st_as_sf(fdata,coords=c("lon","lat"),crs=4326)
  
  # Define thematic map
  # First - the background satellite image
  # Second - the grey "burnt area" hotspots
  # Last - the coloured current hotspots
  m = tm_shape(syd) + tm_rgb()+
    tm_shape(fdata) + tm_squares(alpha=0.5,size=0.03,border.lwd=NA,col = "grey40") +
    tm_shape(tdata) + tm_squares(alpha=0.5,col="max",size=.03,border.lwd=NA,breaks=c(0,50,100,200,400,800,Inf),palette=c("brown","darkred","red","orange","yellow","white")) 
  
  # Add background colour, title etc.
  m = m + tm_layout(bg.color = "grey20",legend.show = FALSE,title = as.character(unq_times[tt]),title.color = "white")
  
  # Save the map with a frame number to the output directory
  fn = paste0("output/",formatC(tt,6,flag = "0"),".png")
  tmap_save(m,fn,width=800, units="px",dpi=150,type="cairo-png")
  return(tt)
}


# The main animation loop - run that function in parallel over each timestep
doit = future_map_chr(seq_along(unq_times),proc)

# To turn this all into an animation i use the ffmpeg tool from the command line
# This command works well if you run it in the /output/ directory:

# ffmpeg -framerate 24 -i %07d.png -framerate 24 -c:v libx264 -pix_fmt yuv420p video.mp4

Ta

@ozjimbob how do we determine the fire intensity, from this data. It seems there are a lot of spots that should not be considered to be fires?

@dicook that’s something we’re still trying to work out and validate. Certainly intensity values < 100 can probably be filtered out, and I am also experimenting with a “persistence” factor in using this data (eg. A hotspot has to appear in the same place for 3 x 10 minutes to be regarded as a potential fire).

Ok