ozjimbob · January 7, 2020 22:19 · dicook · Jan 15, 2020 · ozjimbob · Jan 15, 2020
diff --git a/him_anim.r b/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


	# 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