Observations coverage (ISPD4.7)

../../_images/ISPD.png

This is a static, poster version of this video. The idea is to show the change in observations coverage with time.

The vertical axis is latitude (90S to 90N at 1 degree resolution), the horizontal axis is time (Jan 1836 to Dec 2015, at 1 month resolution), each vertical line is at a randomly sampled 1-degree of longitude. The points are bright yellow if the corresponding lat:lon:month contains more observations than analysis periods (more than one observation/6 hours), pale yellow if it has some observations but fewer than one per analysis period, and dark grey if it contains no observations.

Code to make the poster

Script to download the observations:

#!/usr/bin/env python

# Retrieve all the v3 observations files from NERSC

import datetime
import IRData.twcr as twcr

for year in range(1836,2016):
    twcr.fetch('observations',datetime.datetime(year,1,1),version='3')

Script to extract the observations at each sampled degree of lon:month. This is very slow to run, and should be parallelised, but I don’t expect to run it again anytime soon so I haven’t bothered.

#!/usr/bin/env python

# For each month in 1836-2015, pick a random 1-degree longitude band
#  and get all the observations in that band for that month.

import datetime
import IRData.twcr as twcr
import pickle
import pandas
import numpy
import os
#import os.path

def next_month(dt0):
    dt1 = dt0.replace(day=1)
    dt2 = dt1 + datetime.timedelta(days=32)
    dt3 = dt2.replace(day=1)
    return dt3

obs=[]
for year in range(1836,2015):
    for month in range(1,13):
        fname= "%s/ISPD_poster/%04d/%02d.pkl" % (os.getenv('SCRATCH'),
                                                 year,month)
        if os.path.exists(fname): continue
        sdate = datetime.datetime(year,month,1,0,0)
        edate = next_month(sdate)
        m_ob = twcr.load_observations(sdate,edate,version='3')
        rand_l = numpy.random.randint(0,360)
        longs=numpy.array(m_ob['Longitude'])
        lats=numpy.array(m_ob['Latitude'])
        w=((numpy.isfinite(longs)) & (numpy.isfinite(lats)) &
                    (longs>=rand_l) & (longs <(rand_l+1)) &
                    (lats>=-90) & (lats <=90))
        s_ob = m_ob[w]
        fname= "%s/ISPD_poster/%04d/%02d.pkl" % (os.getenv('SCRATCH'),
                                                 year,month)
        if not os.path.isdir(os.path.dirname(fname)):
            os.makedirs(os.path.dirname(fname))
        pickle.dump( s_ob, open( fname, "wb" ) )

        

Script to plot the poster:

#!/usr/bin/env python

# Make a poster showing ISPD4.7 obs coverage.

import os
import numpy
import datetime

import iris
import matplotlib
from matplotlib.backends.backend_agg import FigureCanvasAgg as FigureCanvas
from matplotlib.figure import Figure
from matplotlib.patches import Rectangle

import pickle

def next_month(dt0):
    dt1 = dt0.replace(day=1)
    dt2 = dt1 + datetime.timedelta(days=32)
    dt3 = dt2.replace(day=1)
    return dt3

start_year=1836
end_year = 2015

# Plot the images
fig=Figure(figsize=(72,9),              # Width, Height (inches)
           dpi=300,
           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)
# Attach a canvas
canvas=FigureCanvas(fig)
matplotlib.rc('image',aspect='auto')


#ax.add_patch(Rectangle((xmin,ymin),width,height,
#                        facecolor='#708090',fill=True,zorder=1))

xmin=datetime.datetime(start_year,1,1,0)-datetime.timedelta(days=30)
xmax=datetime.datetime(end_year,12,13,23)+datetime.timedelta(days=120)
width=xmax-xmin
ymin=-90
ymax=90
height=ymax-ymin

# For each month, load and plot the observations
def y_to_j(y):
    return numpy.minimum(height-1,numpy.maximum(0, 
            numpy.floor((y-ymin)/(ymax-ymin)*(height-1)))).astype(int)

n_obs=numpy.zeros([(end_year+1-start_year)*12,180])
px=[]
for year in range(start_year,end_year+1):
    for month in range(1,13):
        m_count = (year-start_year)*12 + month-1
        px.append(datetime.datetime(year,month,15))
        fname= "%s/ISPD_poster/%04d/%02d.pkl" % (os.getenv('SCRATCH'),
                                                 year,month)
        s_ob = pickle.load( open( fname, "rb" ) )
        sdate = datetime.datetime(year,month,1,0,0)
        edate = next_month(sdate)
        n_steps = int((edate-sdate).total_seconds()/(3600*6))
        lats=numpy.array(s_ob['Latitude'])
        lat_i=y_to_j(lats)
        for i in range(len(lat_i)):
            n_obs[m_count,lat_i[i]] += 1
        n_obs[m_count,:] /= n_steps

s=n_obs.shape
ax2 = fig.add_axes([0,0,1,1],facecolor='green')
ax2.set_axis_off() # Don't want surrounding x and y axis
nd2=numpy.random.rand(s[1],s[0])
clrs=[]
for shade in numpy.linspace(.22+.01,.16+.01):
    clrs.append((shade,shade,shade,1))
y = numpy.linspace(0,1,s[1])
x = numpy.linspace(0,1,s[0])
img = ax2.pcolormesh(x,y,nd2,
                        cmap=matplotlib.colors.ListedColormap(clrs),
                        alpha=1.0,
                        shading='gouraud',
                        zorder=10)

ax = fig.add_axes([0,0,1,1],facecolor='black')
ax.set_axis_off() # Don't want surrounding x and y axis

ax.set_xlim(xmin,xmax)
ax.set_ylim(ymin,ymax)
ax.set_aspect('auto')
# Plot the observations locations
for i in range(s[0]):
    for j in range(s[1]):
        if n_obs[i,j]==0: continue
        scale=min(1.0,max(0.25,numpy.sqrt(n_obs[i,j])))
        s_colour=(.19*(1-scale)+1*scale,
                  .19*(1-scale)+0.84*scale,
                  .19*(1-scale)+0*scale,
                  1)
        ax.add_patch(Rectangle((px[i],j/(s[1]+2)*height+ymin+1),
                               datetime.timedelta(days=31),
                               (height/180),
                               facecolor=s_colour,
                               edgecolor=s_colour,
                               linewidth=0.001,
                               fill=True,
                               zorder=100,
                               alpha=1))
        

fig.savefig('ISPD.pdf')