Global observations coverage video¶
A sample of the observations available for reconstructing past climate. These are pressure observations in the International Surface Pressure Databank (version 3.2.9).
Each observation is shown on-screen for 6-hours either side of the time of observation, and the video shows every observation available, sampling one month every 5 years over the period 1851-2000. The idea is to give an indication of the observations coverage and how it changes with time: In times and places with a fixed red dot we have a continuous series of observations with good diurnal sampling (at least 4 times a day); flashing red dots indicate poor diurnal sampling (usually only once or twice a day), and in too many places there are no observations at all.
The video runs at 12 hours to the second, so to show the full 163 years covered by this dataset would take nearly three days. Instead it shows a sample: all the data for one month, and then a 5-year step to the next month shown. Even so it’s probably best to view selected periods rather than watching from start to end.
Code to make the figure¶
Download the data required:
# Retrieve ISPD3.2.9 observations from NERSC
import datetime
import IRData.twcr as twcr
for year in range(1851,2001):
dte=datetime.datetime(year,1,1)
twcr.fetch_observations(dte,version='2c')
Script to make an individual frame - takes year, month, day, hour, and minute as command-line options:
#!/usr/bin/env python
# Plot the observations in ISPD3.2.9 that are active in 20CR
# for a given time.
import IRData.twcr as twcr
import datetime
import numpy
import pandas
import Meteorographica as mg
import iris
import os
import matplotlib
from matplotlib.backends.backend_agg import FigureCanvasAgg as FigureCanvas
from matplotlib.figure import Figure
import cartopy
import cartopy.crs as ccrs
# Get the datetime to plot from commandline arguments
import argparse
parser = argparse.ArgumentParser()
parser.add_argument("--year", help="Year",
type=int,required=True)
parser.add_argument("--month", help="Integer month",
type=int,required=True)
parser.add_argument("--day", help="Day of month",
type=int,required=True)
parser.add_argument("--hour", help="Hour of day (0 to 23)",
type=int,required=True)
parser.add_argument("--minute", help="Minute of hour (0 to 59)",
type=int,required=True)
parser.add_argument("--radius", help="Marker size (degrees)",
type=float,default=0.75)
parser.add_argument("--opdir", help="Directory for output files",
default="%s/images/transcription_figures/ispd" % \
os.getenv('SCRATCH'),
type=str,required=False)
args = parser.parse_args()
if not os.path.isdir(args.opdir):
os.makedirs(args.opdir)
dte=datetime.datetime(args.year,args.month,args.day,
args.hour,args.minute)
# Define the figure (page size, background color, resolution, ...
aspect=16/9.0
fig=Figure(figsize=(10.8*aspect,10.8), # Width, Height (inches)
dpi=100,
facecolor=(0.88,0.88,0.88,1),
edgecolor=None,
linewidth=0.0,
frameon=False, # Don't draw a frame
subplotpars=None,
tight_layout=None)
# Attach a canvas
canvas=FigureCanvas(fig)
# All mg plots use Rotated Pole, in this case just use the standard
# pole location.
projection=ccrs.RotatedPole(pole_longitude=180.0, pole_latitude=90.0)
# Define an axes to contain the plot. In this case our axes covers
# the whole figure
ax = fig.add_axes([0,0,1,1],projection=projection)
ax.set_axis_off() # Don't want surrounding x and y axis
# Set the axes background colour
ax.background_patch.set_facecolor((0.88,0.88,0.88,1))
# Lat and lon range (in rotated-pole coordinates) for plot
extent=[-180.0,180.0,-90.0,90.0]
ax.set_extent(extent, crs=projection)
# Lat:Lon aspect does not match the plot aspect, ignore this and
# fill the figure with the plot.
matplotlib.rc('image',aspect='auto')
# Draw a lat:lon grid
mg.background.add_grid(ax,
sep_major=5,
sep_minor=2.5,
color=(0,0.3,0,0.2))
# Add the land
land_img=ax.background_img(name='GreyT', resolution='low')
# Load the observations within 6 hours
obs=twcr.load_observations(dte-datetime.timedelta(hours=12),
dte+datetime.timedelta(hours=12),version='2c')
# Time difference in s
dtm=pandas.to_datetime(obs.UID.str.slice(0,10),format="%Y%m%d%H")
dts=numpy.divide(dte-dtm,numpy.timedelta64(1, 's'))
# weight the obs by distance from current time
weight=numpy.abs(dts)/(3600.0*6)
obs['weight']=weight
obs=obs[weight<1]
obs['weight']=1-obs['weight']
# Plot the observations
mg.observations.plot(ax,obs,radius=args.radius,
facecolor='red')
# Add a label showing the date
mg.utils.plot_label(ax,
('%04d-%02d-%02d:%02d' %
(args.year,args.month,args.day,args.hour)),
facecolor=fig.get_facecolor(),
x_fraction=0.99,
y_fraction=0.98,
horizontalalignment='right',
verticalalignment='top',
fontsize=14)
# Render the figure as a png
fig.savefig('%s/obs_%04d%02d%02d%02d%02d.png' %
(args.opdir,args.year,args.month,args.day,
args.hour,args.minute))
To make the video, it is necessary to run the script above hundreds of times - giving an image for every 30-minute period in the sampled months. 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 time
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)
# Plotted size of observations might vary with date
def obs_size(dte):
return 0.75
# Timestep might vary with date as well
# Also need to put in some big skips, to keep
# total length somewhere reasonable.
def next_timestep(dte):
old=dte
new=dte+datetime.timedelta(minutes=30)
if(old.month!=new.month):
new=datetime.datetime(new.year+5,
new.month,new.day,
new.hour,new.minute)
return new
start_day=datetime.datetime(1851, 1, 1,13)
end_day =datetime.datetime(2000,12,31,11)
# Function to check if the job is already done for this timepoint
def is_done(dte):
op_file_name=("%s/images/transcription_figures/ispd/" +
"obs_%04d%02d%02d%02d%02d.png") % (
os.getenv('SCRATCH'),
dte.year,dte.month,dte.day,
dte.hour,dte.minute)
if os.path.isfile(op_file_name):
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:
if is_done(current_day):
current_day=next_timestep(current_day)
continue
f=open("multirun.slm","w+")
f.write('#!/bin/ksh -l\n')
f.write('#SBATCH --output=%s/TR_obs_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=10\n')
for i in (0,1,2,3):
while is_done(current_day):
current_day=next_timestep(current_day)
if current_day>=end_day: break
cmd=("./plot_obs_fortime.py --year=%d --month=%d" +
" --day=%d --hour=%d --minute=%d & \n") % (
current_day.year,current_day.month,
current_day.day,current_day.hour,
current_day.minute)
f.write(cmd)
current_day=next_timestep(current_day)
f.write('wait')
f.close()
rc=subprocess.call('sbatch multirun.slm',shell=True)
os.unlink('multirun.slm')
#sys.exit(0)
max_new_jobs=max_new_jobs-1
time.sleep(60)
To turn the thousands of images into a movie, use ffmpeg
ffmpeg -r 24 -pattern_type glob -i ispd/\*.png \
-c:v libx264 -threads 16 -preset slow -tune animation \
-profile:v high -level 4.2 -pix_fmt yuv420p -crf 25 \
-c:a copy ispd.mp4