Observations coverage video¶
In each 1x1 degree grid-cell, a bright yellow circle is shown if at least one pressure observation is available every 6 hours (one in each assimilation run). Paler yellow circles indicate observations in some, but not all assimilation periods (partial coverage). Each frame covers a 10-day period.
Code to make the figure¶
Script to make an individual frame - takes year, month, day, and hour as command-line options:
#!/usr/bin/env python
# Observations coverage in 20CRv3
import os
import datetime
import iris
import numpy
import pandas
import matplotlib
from matplotlib.backends.backend_agg import FigureCanvasAgg as FigureCanvas
from matplotlib.figure import Figure
from matplotlib.patches import Rectangle
from matplotlib.lines import Line2D
# Fix dask SPICE bug
import dask
dask.config.set(scheduler='single-threaded')
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="Time of day (0 to 23.99)",
default=12,
type=float,required=False)
parser.add_argument("--opdir", help="Directory for output files",
default="%s/images/20CRv3_observations" % \
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,
int(args.hour),int(args.hour%1*60))
def load_observations_1file(of_name):
compression='infer'
if of_name[-2:]=='gz': compression='gzip'
o=pandas.read_fwf(of_name,
colspecs=[(0,19),
(20,50),
(52,64),
(66,68),
(69,72),
(74,80),
(81,87),
(88,95),
(97,102),
(103,110),
(111,112),
(113,123),
(124,134),
(135,145),
(146,156),
(157,167),
(168,175),
(176,183),
(184,191),
(192,200),
(201,205),
(206,213),
(214,221),
(222,223)],
header=None,
encoding="ISO-8859-1",
compression=compression,
usecols=[0,1,2,3,4,5,6],
names=['UID',
'Name',
'ID',
'Type',
'NCEP.Type',
'Longitude',
'Latitude',
'Observed',
'Time.offset',
'Observed.2',
'Skipped',
'Bias.correction',
'Obfit.prior',
'Obfit.post',
'Obsprd.prior',
'Obsprd.post',
'Oberrvar.orig.out',
'Oberrvar.out',
'Oberrvar.use',
'Paoverpb.save',
'Prob.gross.error',
'Localization.length.scale',
'Lnsigl',
'QC.failure.flag'],
converters={'UID': str,
'Name': str,
'ID': str,
'Type': str,
'NCEP.Type': str,
'Longitude': float,
'Latitude': float,
'Observed': float,
'Time.offset': float,
'Observed.2': float,
'Skipped': int,
'Bias.correction': float,
'Obfit.prior': float,
'Obfit.post': float,
'Obsprd.prior': float,
'Obsprd.post': float,
'Oberrvar.orig.out': float,
'Oberrvar.out': float,
'Oberrvar.use': float,
'Paoverpb.save': float,
'Prob.gross.error': float,
'Localization.length.scale': float,
'Lnsigl': float,
'QC.failure.flag': int},
na_values=['NA','*','***','*****','*******','**********',
'-99','9999','-999','9999.99','10000.0',
'-9.99',
'999999999999','9'],
comment=None)
return(o)
# Load the obs for +-15 days around given datetime
# Get the fraction of assimilation points with at least
# 1 ob for each 1x1 degree grid-cell.
n_fields=0
width=360
height=180
xmin=-180
xmax=180
ymin=-90
ymax=90
n_obs=numpy.zeros([width,height])
def x_to_i(x):
return numpy.minimum(width-1,numpy.maximum(0,
numpy.floor((x-xmin)/(xmax-xmin)*(width-1)))).astype(int)
def y_to_j(y):
return numpy.minimum(height-1,numpy.maximum(0,
numpy.floor((y-ymin)/(ymax-ymin)*(height-1)))).astype(int)
def i_to_x(i):
return xmin + ((i+1)/width) * (xmax-xmin)
def j_to_y(j):
return ymin + ((j+1)/height) * (ymax-ymin)
ct=dte-datetime.timedelta(days=15)
et=dte+datetime.timedelta(days=15)
while ct<et:
if ct.hour%6==0:
ofile="%s/20CR/version_3/v3_observations/%04d/%04d%02d%02d%02d_psobs_posterior.txt" % (
os.getenv('SCRATCH'),ct.year,ct.year,ct.month,ct.day,ct.hour)
if not os.path.isfile(ofile):
ofile+='.gz'
obs=load_observations_1file(ofile)
longs=numpy.array(obs['Longitude'])
lats=numpy.array(obs['Latitude'])
w=((numpy.isfinite(longs)) & (numpy.isfinite(lats)) &
(longs>=xmin) & (longs <=360) &
(lats>=ymin) & (lats <=ymax))
longs=longs[w]
lats=lats[w]
longs[longs>180] -= 360
lon_i=x_to_i(longs)
lat_i=y_to_j(lats)
n_add=numpy.zeros([width,height])
for i in range(len(lon_i)):
n_add[lon_i[i],lat_i[i]] = 1
n_obs += n_add
n_fields += 1
ct += datetime.timedelta(hours=1)
n_obs /= n_fields
# Load a land mask
mask = iris.load_cube("%s/fixed_fields/land_mask/opfc_global_2019.nc" % os.getenv('DATADIR'))
# Define the figure (page size, background color, resolution, ...
fig=Figure(figsize=(19.2,10.8), # Width, Height (inches)
dpi=100,
facecolor=(0.5,0.5,0.5,1),
edgecolor=None,
linewidth=0.0,
frameon=False, # Don't draw a frame
subplotpars=None,
tight_layout=None)
fig.set_frameon(False)
# Attach a canvas
canvas=FigureCanvas(fig)
# Projection for plotting
cs=iris.coord_systems.RotatedGeogCS(90,180,0)
def plot_cube(resolution):
lat_values=numpy.arange(ymin,ymax+resolution,resolution)
latitude = iris.coords.DimCoord(lat_values,
standard_name='latitude',
units='degrees_north',
coord_system=cs)
lon_values=numpy.arange(xmin,xmax+resolution,resolution)
longitude = iris.coords.DimCoord(lon_values,
standard_name='longitude',
units='degrees_east',
coord_system=cs)
dummy_data = numpy.zeros((len(lat_values), len(lon_values)))
plot_cube = iris.cube.Cube(dummy_data,
dim_coords_and_dims=[(latitude, 0),
(longitude, 1)])
return plot_cube
# Define an axes to contain the plot. In this case our axes covers
# the whole figure
ax = fig.add_axes([0,0,1,1])
ax.set_axis_off() # Don't want surrounding x and y axis
# Lat and lon range (in rotated-pole coordinates) for plot
ax.set_xlim(xmin,xmax)
ax.set_ylim(ymin,ymax)
ax.set_aspect('auto')
# Background
ax.add_patch(Rectangle((xmin,ymin),width,height,facecolor=(0.5,0.5,0.5,1),fill=True,zorder=1))
# Draw lines of latitude and longitude
for lat in range(-90,95,5):
lwd=0.25
x=[]
y=[]
for lon in range(-180,181,1):
rp=iris.analysis.cartography.rotate_pole(numpy.array(lon),
numpy.array(lat),
180,
90)
nx=rp[0]
if nx>180: nx -= 360
ny=rp[1]
if(len(x)==0 or (abs(nx-x[-1])<10 and abs(ny-y[-1])<10)):
x.append(nx)
y.append(ny)
else:
ax.add_line(Line2D(x, y, linewidth=lwd, color=(0.4,0.4,0.4,1),
zorder=10))
x=[]
y=[]
if(len(x)>1):
ax.add_line(Line2D(x, y, linewidth=lwd, color=(0.4,0.4,0.4,1),
zorder=10))
for lon in range(-180,185,5):
lwd=0.25
x=[]
y=[]
for lat in range(-90,90,1):
rp=iris.analysis.cartography.rotate_pole(numpy.array(lon),
numpy.array(lat),
180,
90)
nx=rp[0]
if nx>180: nx -= 360
ny=rp[1]
if(len(x)==0 or (abs(nx-x[-1])<10 and abs(ny-y[-1])<10)):
x.append(nx)
y.append(ny)
else:
ax.add_line(Line2D(x, y, linewidth=lwd, color=(0.4,0.4,0.4,1),
zorder=10))
x=[]
y=[]
if(len(x)>1):
ax.add_line(Line2D(x, y, linewidth=lwd, color=(0.4,0.4,0.4,1),
zorder=10))
# Plot the land mask
mask_pc=plot_cube(0.05)
mask = mask.regrid(mask_pc,iris.analysis.Linear())
lats = mask.coord('latitude').points
lons = mask.coord('longitude').points
mask_img = ax.pcolorfast(lons, lats, mask.data,
cmap=matplotlib.colors.ListedColormap(
((0,0,0,0),
(0,0,0,1))),
vmin=0,
vmax=1,
alpha=1.0,
zorder=20)
# Plot the observations locations
for i in range(width):
for j in range(height):
if n_obs[i,j]==0: continue
rp=iris.analysis.cartography.rotate_pole(numpy.array(i_to_x(i)),
numpy.array(j_to_y(j)),
180,
90)
nlon=rp[0][0]
nlat=rp[1][0]
ax.add_patch(matplotlib.patches.Circle((nlon,nlat),
radius=0.49,
facecolor='yellow',
edgecolor='yellow',
linewidth=0.1,
alpha=max(0.15,n_obs[i,j]),
zorder=180))
# Label with the date
ax.text(xmax-width*0.009,
ymax-height*0.016,
"%04d-%02d" % (args.year,args.month),
horizontalalignment='right',
verticalalignment='top',
color='black',
bbox=dict(facecolor=(0.8,0.8,0.8,0.5),
edgecolor='black',
boxstyle='round',
pad=0.5),
size=14,
clip_on=True,
zorder=500)
# Render the figure as a png
fig.savefig('%s/%04d%02d%02d.png' % (args.opdir,args.year,
args.month,args.day))
To make the video, it is necessary to run the script above hundreds of times - giving an image for every 10-day period. This script makes the list of commands needed to make all the images, which can be run in parallel.
#!/usr/bin/env python
# Make all the individual frames for a movie
import os
import subprocess
import datetime
# Where to put the output files
opdir="%s/slurm_output" % os.getenv('SCRATCH')
if not os.path.isdir(opdir):
os.makedirs(opdir)
# Function to check if the job is already done for this timepoint
def is_done(year,month,day):
op_file_name=("%s/images/20CRv3_observations/" +
"%04d%02d%02d.png") % (
os.getenv('SCRATCH'),
year,month,day)
if os.path.isfile(op_file_name):
return True
return False
f=open("run.txt","w+")
start_day=datetime.datetime(1851, 1, 1, 0)
end_day =datetime.datetime(2015, 12, 31, 23)
current_day=start_day
while current_day<=end_day:
if is_done(current_day.year,current_day.month,
current_day.day):
current_day=current_day+datetime.timedelta(days=10)
continue
cmd=("./plot_obs_coverage.py --year=%d --month=%d " +
"--day=%d "+
"\n") % (
current_day.year,current_day.month,
current_day.day)
f.write(cmd)
current_day=current_day+datetime.timedelta(days=10)
f.close()
To turn the thousands of images into a movie, use ffmpeg
ffmpeg -r 24 -pattern_type glob -i 20CRv3_observations/\*.png \
-c:v libx264 -threads 16 -preset veryslow -tune film \
-profile:v high -level 4.2 -pix_fmt yuv420p \
-b:v 5M -maxrate 5M -bufsize 20M \
-c:a copy 20CRv3_observations.mp4
