cd /Users/laban/devel/git/private/forecast-weather


from scipy.io import netcdf
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
import numpy as np
def head(x): return x[:10]; end

data = netcdf.netcdf_file("./land.netcdf")

data.variables

lat = data.variables['latitude']
lon = data.variables['longitude']
surface = data.variables['LAND_surface']
p_time = data.variables['time']

import time
print "Data for:", time.ctime(p_time.data[0])

print surface.data[0]
surface.data[0].shape
surface_bitmap = surface.data[0] * 255

fig = plt.figure(figsize=(12, 6))
plt.imshow(surface_bitmap, origin='lower')

fig = plt.figure(figsize=(12, 6))
plt.contour(surface_bitmap, origin='lower', colors="blue", aspect='auto', linewidths=0.1)

temp_data = netcdf.netcdf_file("temp2.netcdf").variables['TMP_2maboveground']
print temp_data.shape
temp_2m_above_ground_deg_C = temp_data[0] - 273.15


fig = plt.figure(figsize=(12, 6))
plt.contour(surface_bitmap, origin='lower', colors="black", aspect='auto', linewidths=0.1)
plt.imshow(temp_2m_above_ground_deg_C, origin='lower')
colorbar()

africa_coords = [slice(250,700), slice(0,250)]
temp_2m_above_ground_deg_C[africa_coords].shape
africa_temp = temp_2m_above_ground_deg_C[africa_coords]

# Sahara's scorching
fig = plt.figure(figsize=(12, 6))
plt.imshow(africa_temp, origin='lower')
colorbar()

veg_data = netcdf.netcdf_file("./vegetation.netcdf")
vegetation = veg_data.variables['VEG_surface'][0] 
world_vegetation = vegetation.copy()

# Truncate values that are excessively large
world_vegetation[world_vegetation > 1e3] = 0




# Central Africa is the only significant bit of vegetation left in Africa :(
fig = plt.figure(figsize=(12, 6))
world_vegetation[africa_coords].shape
africa_vegetation = world_vegetation[africa_coords]
plt.imshow(africa_vegetation, origin='lower', cmap=get_cmap('BuGn'))
colorbar()

# Amazon's still green. Yay!
fig = plt.figure(figsize=(12, 6))
plt.imshow(world_vegetation, origin='lower', cmap=get_cmap('BuGn'))
colorbar()

# Cloud cover

cloud_data = netcdf.netcdf_file("./cloud.netcdf")
clouds = cloud_data.variables['TCDC_lowcloudlayer'][0] 


fig = plt.figure(figsize=(12, 6))
ax = fig.add_subplot(1,1,1)
plt.contour(surface_bitmap, origin='lower', colors="black", aspect='auto', linewidths=0.4)
im = plt.imshow(clouds, origin='lower', cmap=get_cmap('BuGn'))

# Let's put in some lat-long on the cloudy data

# X ticks
xticks = np.linspace(0,lon.data.shape[0], 9)
xlabels = np.linspace(0,315, 8)
ax.set_xticks(xticks)
ax.set_xticklabels(xlabels) 

# Y ticks
yticks = np.linspace(0,lat.data.shape[0], 5)
ylabels = np.linspace(-90,90, 5)
ax.set_yticks(yticks)
ax.set_yticklabels(ylabels) 




africa_clouds = clouds[africa_coords]
africa_contour = surface_bitmap[africa_coords]
plt.contour(africa_contour, origin='lower', colors="blue", aspect='auto', linewidths=0.1)
plt.imshow(africa_clouds, origin='lower', cmap=get_cmap('YlGn'))
colorbar()

# Horn of Africa
poi = [slice(430,500), slice(160,260)]
poi_contour = surface_bitmap[poi]

# plot it
# 
fig = plt.figure(figsize=(12, 6))

ax = plt.subplot(1,3,1)
plt.title("% cloud cover")
poi_clouds = clouds[poi]
plt.contour(poi_contour, origin='lower', colors="blue", aspect='auto', linewidths=0.1)
plt.imshow(poi_clouds, origin='lower', cmap=get_cmap('YlGn'))
colorbar()

ax = plt.subplot(1,3,2)
plt.title("Temp in C")
temp_poi = temp_2m_above_ground_deg_C[poi]
plt.contour(poi_contour, origin='lower', colors="blue", aspect='auto', linewidths=0.1)
plt.imshow(temp_poi, origin='lower')
colorbar()

ax = plt.subplot(1,3,3)
plt.title("% vegetation cover")
poi_vegetation = world_vegetation[poi]
plt.contour(poi_contour, origin='lower', colors="blue", aspect='auto', linewidths=0.1)
plt.imshow(poi_vegetation, origin='lower', cmap=get_cmap('BuGn'))
colorbar()

# Colour map from www.scipy.org/Cookbook/Matplotlib/Show_colormaps
from pylab import *
from numpy import outer
rc('text', usetex=False)
a=outer(arange(0,1,0.01),ones(10))
figure(figsize=(10,5))
subplots_adjust(top=0.8,bottom=0.05,left=0.01,right=0.99)
maps=[m for m in cm.datad if not m.endswith("_r")]
maps.sort()
l=len(maps)+1
for i, m in enumerate(maps):
    subplot(1,l,i+1)
    axis("off")
    imshow(a,aspect='auto',cmap=get_cmap(m))
    title(m,rotation=90,fontsize=10)
#savefig("colormaps.png",dpi=100,facecolor='gray')
?ax

# And that's it folks!