Analysis of insolation and air temperature (thread: http://www.reddit.com/r/meteorology/comments/2t01gt/question_whats_going_with_heat_insolation/)

analysis.py

import pandas as pd
import numpy as np
import re
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
 
def load_grid(path):
    f = open(path, 'r')
    
    headers = None
    for line in f:
        if line.startswith('Lat Lon'):
            headers = re.split('\s+', line.strip())
            break
 
    records = []
    for line in f:
        records.append([float(x) for x in re.split('\s+', line.strip())])
 
    df = pd.DataFrame(np.array(records), columns=headers)
    assert len(df) == len(set(zip(df['Lat'], df['Lon'])))
 
    return df
 
 
# https://eosweb.larc.nasa.gov/cgi-bin/sse/global.cgi?email=skip@larc.nasa.gov
 
# Insolation Incident On A Horizontal Surface (kWh/m^2/day)
horizontal_radiation = load_grid('Horizontal_Radiation.txt')
# Air Temperature at 10 m Above The Surface Of The Earth (deg C)
air_temp = load_grid('AirTemp_22yr_T10M.txt')
 
insolation_efficiency = []
 
assert np.all(horizontal_radiation['Lon'] == air_temp['Lon'])
assert np.all(horizontal_radiation['Lat'] == air_temp['Lat'])
 
for i in xrange(len(horizontal_radiation)):
    if horizontal_radiation.loc[i,'Lat'] >= 0:
        stat = (air_temp.iloc[i,8]+273.0) / horizontal_radiation.iloc[i,8]
    else:
        stat = (air_temp.iloc[i,2]+273.0) / horizontal_radiation.iloc[i,2]
 
    insolation_efficiency.append((
        horizontal_radiation.loc[i,'Lat'],
        horizontal_radiation.loc[i,'Lon'],
        stat))
 
insolation_efficiency = pd.DataFrame(np.array(insolation_efficiency), columns=['Lat', 'Lon', 'ieff'])
 
lats = np.arange(-90,90,1)
lons = np.arange(-180,180,1)
lons, lats = np.meshgrid(lons,lats)
 
fig = plt.figure()
ax = fig.add_axes([0.05,0.05,0.9,0.9])
 
m = Basemap(projection='kav7',lon_0=0,resolution='l')
 
m.drawcoastlines(linewidth=0.25)
m.drawcountries(linewidth=0.25)
m.drawmapboundary(fill_color='#A0D9F1')
 
data = insolation_efficiency['ieff'].values
data = data.reshape((180, 360))
 
im1 = m.pcolormesh(
    lons, lats,
    data,
    shading='flat',
    cmap=plt.cm.jet,
    latlon=True)
 
cb = m.colorbar(im1,"bottom", size="5%", pad="2%")
 
# draw parallels and meridians, but don't bother labelling them.
m.drawparallels(np.arange(-90.,99.,30.))
m.drawmeridians(np.arange(-180.,180.,60.))
 
ax.set_title('Air Temperature @10m ($K$) / Horizontal Insolation ($kWh/m^2/day$) in Summer Month')
plt.show()