20CR3 compared with new observations in January 1872 (video)ΒΆ
Left panel is a spaghetti-contour plot of reanalysis mean-sea-level pressure. The small yellow dots mark stations assimilated by the reanalysis; the larger grey dots mark the new stations used for validation.
Right panel shows the observed MSLP (red line) and the reanalysis ensemble values (blue dots) at the station location, for each new station.
Download the data (uses this script)
#!/bin/sh
../../../scripts/get_data.py --reanalysis=20cr3 --start="1872010118" --end="1872013118"
To make the video, it is necessary to run this script hundreds of times - giving an image for every 15-minute period. The best way to do this is system dependent - the script below does it on the Met Office SPICE cluster - it will need modification to run on any other system. (Could do this on a single PC, but it will take many hours).
#!/usr/bin/env python
# Make all the individual frames for a movie
# run the jobs on SPICE.
import os
import sys
import subprocess
import datetime
max_jobs_in_queue=500
# Where to put the output files
opdir="%s/slurm_output" % os.getenv('SCRATCH')
if not os.path.isdir(opdir):
os.makedirs(opdir)
start_day=datetime.datetime(1872, 1, 1, 13)
end_day =datetime.datetime(1872, 1, 31, 11)
# Output directory
opd="%s/images/DWR/vcs_20CR3_1872_scatter+contour" % os.getenv('SCRATCH')
if not os.path.isdir(opd):
os.makedirs(opd)
# Get file name from date
def opfile(dte):
return "%s/Scatter+contour_%04d%02d%02d%02d%02d.png" % (
opd,dte.year,dte.month,dte.day,dte.hour,dte.minute)
# Function to check if the job is already done for this timepoint
def is_done(dte):
if os.path.isfile(opfile(dte)):
return True
return False
current_day=start_day
while current_day<=end_day:
queued_jobs=subprocess.check_output('squeue --user hadpb',
shell=True).count('\n')
max_new_jobs=max_jobs_in_queue-queued_jobs
while max_new_jobs>0 and current_day<=end_day:
f=open("multirun.slm","w+")
f.write('#!/bin/ksh -l\n')
f.write(("#SBATCH --output=%s/sc_frame_"+
"%04d%02d%02d%02d.out\n") %
(opdir,
current_day.year,current_day.month,
current_day.day,current_day.hour))
f.write('#SBATCH --qos=normal\n')
f.write('#SBATCH --ntasks=4\n')
f.write('#SBATCH --ntasks-per-core=1\n')
f.write('#SBATCH --mem=40000\n')
f.write('#SBATCH --time=5\n')
count=0
for minute in (0,15,30,45):
if is_done(datetime.datetime(current_day.year,
current_day.month,
current_day.day,
current_day.hour,
minute)):
continue
cmd=("../../../../scripts/scatter+contour.py "+
" --year=%d --month=%d "+
"--day=%d --hour=%f "+
"--opdir=%s "+
"--reanalysis=20cr3 "+
"--video "+
"--skip=ABERDEEN "+
"--skip=LISBON "+
"--skip=ROCHEFORT "+
"--skip=BASRA "+
"--skip=DIYARBAKIR "+
"--skip=BAGHDAD &\n") % (
current_day.year,current_day.month,
current_day.day,current_day.hour+minute/60.0,
opd)
f.write(cmd)
count=count+1
f.write('wait\n')
f.close()
current_day=current_day+datetime.timedelta(hours=1)
if count>0:
max_new_jobs=max_new_jobs-1
rc=subprocess.call('sbatch multirun.slm',shell=True)
os.unlink('multirun.slm')
To turn the thousands of images into a movie, use ffmpeg
ffmpeg -r 24 -pattern_type glob -i vcs_20CR3_1872_scatter+contour/\*.png \
-c:v libx264 -threads 16 -preset slow -tune animation \
-profile:v high -level 4.2 -pix_fmt yuv420p -crf 25 \
-c:a copy vcs_20CR3_1872_scatter+contour.mp4
