# Why does the sky look gray through this window?¶

This was prepared as part of my answer to this physics.SE question

in which the user shares two pictures taken out of the same window seconds apart.

Here we will try to validate the claim by Luboš Motl in this answer that it just appears grey because of the additional scattered light from the sun hitting the windows.

In [165]:
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
from skimage.color import convert_colorspace
from mpl_toolkits.axes_grid1 import make_axes_locatable
from mpl_toolkits.axes_grid1 import ImageGrid
%matplotlib inline

In [166]:
bluesky = imread("http://i.stack.imgur.com/fbIyO.jpg")


Next we will convert the images from the fairly useless RGB colorspace to the much more useful XYZ colorspace which is built on a model of the actual receptors in the human eye. The $Y$ coordinate corresponds to the percieved luminance of the image (i.e. the average human response across the visible spectrum) and the $Z$ coordinate corresponds to our blue receptor response. The $X$ coordinate is set to be orthogonal to the other two. See here the responses from wikipedia:

In [167]:
bluesky_xyz = convert_colorspace(bluesky,"RGB","XYZ")
greysky_xyz = convert_colorspace(greysky,"RGB","XYZ")

print bluesky_xyz.min(), bluesky_xyz.max()
print greysky_xyz.min(), greysky_xyz.max()

4.10787348995e-05 1.088754
0.0 1.088754

In [168]:
#better colorscheme
plt.rc("image", cmap="cubehelix")

In [169]:
fig,axs = plt.subplots(1,2, subplot_kw={"xticks":[], "yticks":[]})
axs[0].imshow(bluesky_xyz[:,:,1],vmin=0,vmax=1);
axs[0].set_title("Blue");
axs[1].imshow(greysky_xyz[:,:,1],vmin=0,vmax=1);
axs[1].set_title("Gray");


This certainly seems to support the notion that the luminosity has increased. Let's take a look at the blue channel and see if anything has changed.

In [171]:
fig,axs = plt.subplots(1,2, subplot_kw={"xticks":[], "yticks":[]})
axs[0].imshow(bluesky_xyz[:,:,2],vmin=0,vmax=1);
axs[0].set_title("Blue");
axs[1].imshow(greysky_xyz[:,:,2],vmin=0,vmax=1);
axs[1].set_title("Gray");


Honestly, these look pretty similar to me. If anything the amount of blue in the "Gray" image has increased as well.

Putting them all together to make a pretty picture.

In [172]:
fig = plt.figure(1, (27,8))
grid = ImageGrid(fig, 111, nrows_ncols=(3,2), axes_pad=0.1,
share_all=True, label_mode="L", cbar_location="bottom", cbar_mode="single" )

grid[0].imshow(bluesky);
grid[1].imshow(greysky);
grid[2].imshow(bluesky_xyz[:,:,1],vmin=0,vmax=1);
grid[3].imshow(greysky_xyz[:,:,1],vmin=0,vmax=1);
grid[4].imshow(bluesky_xyz[:,:,2],vmin=0,vmax=1);
im = grid[5].imshow(greysky_xyz[:,:,2],vmin=0,vmax=1);
grid.cbar_axes[0].colorbar(im)
grid[0].set_title("Blue")
grid[1].set_title("Gray")
grid[0].set_ylabel("original")
grid[2].set_ylabel("Y = luminance")
grid[4].set_ylabel("Z ~ blue")
grid.axes_llc.set_xticks([]);
grid.axes_llc.set_yticks([]);


Let's next look at the histogram of pixel values for the $Y$ and $Z$ coordinates to see if we can tell what is going on.

In [152]:
fig,axs = plt.subplots(2, figsize=(10,3), sharex=True, sharey=True);
bins = np.linspace(0,1,200)
x = bins[1:] - 0.5*(bins[1]-bins[0])
axs[0].plot(x,np.histogram(bluesky_xyz[:,:,1].ravel(),bins,normed=True)[0],c='gray',lw=2, alpha=0.9, label="Y");
axs[0].plot(x,np.histogram(bluesky_xyz[:,:,2].ravel(),bins,normed=True)[0],c='blue',lw=2, alpha=0.7, label="Z");
axs[1].plot(x,np.histogram(greysky_xyz[:,:,1].ravel(),bins,normed=True)[0],c='gray',lw=2, alpha=0.9);
axs[1].plot(x,np.histogram(greysky_xyz[:,:,2].ravel(),bins,normed=True)[0],c='blue',lw=2, alpha=0.7);
plt.ylim((0,5));
plt.xlim((0,1.0));
axs[0].legend();
plt.yticks([]);
axs[0].vlines([0.3,0.6],0,5,linestyle="--",alpha=0.3);
axs[1].vlines([0.35,0.6],0,5,linestyle="--",alpha=0.3);
axs[0].set_ylabel("Blue Sky");
axs[1].set_ylabel("Gray Sky");


Looking at this histograms, we can see clear band of blues with middle type activations, this is probably the sky, let's check

In [173]:
fig,axs = plt.subplots(1, 2, figsize=(10,5), subplot_kw={"xticks":[], "yticks":[]});
skyfilt_blue = (bluesky_xyz[:,:,2] > 0.3)*(bluesky_xyz[:,:,2] < 0.6)
skyfilt_grey = (greysky_xyz[:,:,2] > 0.35)*(greysky_xyz[:,:,2] < 0.6)
axs[0].imshow( skyfilt_blue, cmap=plt.bone());
axs[1].imshow( skyfilt_grey, cmap=plt.bone());

In [177]:
fig,axs = plt.subplots(1, 2, figsize=(10,5), subplot_kw={"xticks":[], "yticks":[]});
blue_emph = bluesky.copy()
blue_emph[~skyfilt_blue] *= 0
grey_emph = greysky.copy()
grey_emph[~skyfilt_grey] *= 0
axs[0].imshow( blue_emph, cmap=plt.bone());
axs[0].set_title("Blue Sky");
axs[1].set_title("Gray Sky");
axs[1].imshow( grey_emph, cmap=plt.bone());


Yup its the sky, let's look at the filtered histograms.

In [175]:
fig,axs = plt.subplots(2, figsize=(10,3), sharex=True, sharey=True);
bins = np.linspace(0,1,200)
x = bins[1:] - 0.5*(bins[1]-bins[0])
axs[0].plot(x,np.histogram(bluesky_xyz[skyfilt_blue][:,1],bins,normed=True)[0],c='gray',lw=2, alpha=0.9, label="Y");
axs[0].plot(x,np.histogram(bluesky_xyz[skyfilt_blue][:,2],bins,normed=True)[0],c='blue',lw=2, alpha=0.7, label="Z");
axs[1].plot(x,np.histogram(greysky_xyz[skyfilt_grey][:,1],bins,normed=True)[0],c='gray',lw=2, alpha=0.9);
axs[1].plot(x,np.histogram(greysky_xyz[skyfilt_grey][:,2],bins,normed=True)[0],c='blue',lw=2, alpha=0.7);
plt.ylim((0,15));
plt.xlim((0.2,0.6));
plt.yticks([]);
axs[0].legend();
axs[0].vlines([0.3,0.6],0,15,linestyle="--",alpha=0.3);
axs[1].vlines([0.35,0.6],0,15,linestyle="--",alpha=0.3);
axs[0].set_ylabel("Blue Sky");
axs[1].set_ylabel("Gray Sky");

In [176]:
fig,axs = plt.subplots(2, figsize=(10,3), sharex=True, sharey=True);
bins = np.arange(256)
x = bins[:-1]
axs[0].plot(x,np.histogram(bluesky[skyfilt_blue][:,0],bins,normed=True)[0],c='red',lw=2, alpha=0.7, label="R");
axs[0].plot(x,np.histogram(bluesky[skyfilt_blue][:,1],bins,normed=True)[0],c='green',lw=2, alpha=0.7, label="G");
axs[0].plot(x,np.histogram(bluesky[skyfilt_blue][:,2],bins,normed=True)[0],c='blue',lw=2, alpha=0.7, label="B");
axs[1].plot(x,np.histogram(greysky[skyfilt_grey][:,0],bins,normed=True)[0],c='red',lw=2, alpha=0.7);
axs[1].plot(x,np.histogram(greysky[skyfilt_grey][:,1],bins,normed=True)[0],c='green',lw=2, alpha=0.7);
axs[1].plot(x,np.histogram(greysky[skyfilt_grey][:,2],bins,normed=True)[0],c='blue',lw=2, alpha=0.7);
# plt.ylim((0,15));
plt.xlim((100,210));
plt.yticks([]);
axs[0].legend();
axs[0].set_ylabel("Blue Sky");
axs[1].set_ylabel("Gray Sky");

In [ ]: