Aerosol attenuation in the atmosphere has a relatively weak spectral variation (weak compared to molecular absorption). However at high solar zenith angles, with the sun towards the horizon, the spectral irradiance on object at sea level differs substantially from the spectral irradiance with the sun at zero zenith angle. The radiance from an object comprises at least three distinctly different components: the (1) reflected direct sunlight, (2) the reflected sky shine, and (3) the reflected earth shine.
The perceived colour of an object depends on the object's inherent colour as well as the irradiance spectrum. The spectral properties of the three irradiance components on an object at sea level vary dramatically with solar zenith angle. When the sun is at zero zenith angle the sunlight is least attenuated (spectral variation is smaller), and the spectral properties of sky shine is similar to the direct sun irradiance. When the sun is near the horizon the direct sun irradiance shifts towards red, while the sky radiance shifts towards blue. Hence, depending on the orientation of the surface, the colour can shift towards red (facing the sun), towards blue (facing away from the sun), or anywhere in between depending on the ratio of direct sun to reflected sky. It is reasonable to refer to a shift in the colour balance of a scene because of the solar zenith angle. This shift is different for different object surface orientations in the scene, depending on the angle between the surface normal vector and the direction to the sun.
Investigate the effect of the slow spectral variation in sunlight on the colour coordinates of a few objects through a selection of atmospheric conditions. Ignore the reflected sky shine and reflected earth shine, consider only direct solar irradiance. The data used in this investigation is in file DP05.tgz available at
https://github.com/NelisW/ComputationalRadiometry/blob/master/DataPackages/DP05.tgz?raw=true. Study the readme.txt file for more information.
Do the analysis for the following samples: Lettuce, Skin Black, Skin White, Dark Green Leaf, Red Paint, Blue Nitrile Glove, and White Paper, being columns 2, 6, 7, 8, 9, 10, and 11 in the samplesVis.txt data file.
Modtran5 was used to calculate the direct solar irradiance at sea level, and the corresponding atmospheric transmittance to the sun. The Modtran atmospheres selected for this investigation are:
The main differences between the above atmospheres are in are in the aerosol models used; the differences in temperature, pressure and relative humidity vertical profiles have little effect on the visual band transmittance (in the infrared bands these differences are not negligible). The solar irradiance was calculated for each of the above scenarios at zenith angles of $\{0, 45, 60, 80, 88\}$, where zero zenith angle is vertically up and 88 degrees is near the horizon.
Load and plot the samples' spectral reflectance. Show all the samples on the same graph. [2]
Load and plot the four atmosphere's spectral transmittance. Show each atmosphere separately in its own plot, but show transmittance for all zenith angles for each atmosphere in the same plot. Study the plots, compare the transmittance values for the different atmospheres and zenith angles and summarise the important similarities and differences. [4]
Load and plot the four atmosphere's spectral solar irradiance. Show each atmosphere separately in its own plot, but show irradiance for all zenith angles for each atmosphere in the same plot. Study the plots, compare the irradiance values for the different atmospheres and zenith angles and summarise the important similarities and differences. [4]
Plot the CIE coordinates of monochromatic colour from 380 nm to 720 nm. It should be a familiar horseshoe shape. Find a colour version of the CIE diagram on the internet and confirm your calculations. [2]
Present a complete mathematical derivation on how the CIE colour coordinates are calculated from spectral radiance. [3]
Provide a description of how the CIE coordinates can be calculated from the RGB coordinates. [3]
Calculate and plot a sequence of graphs (one graph for each atmosphere) showing the CIE colour coordinates for all the samples, showing the change in colour coordinates with different zenith angles. The CIE plots must show the monochromatic horseshoe curve. The idea is to demonstrate the change in colour with varying zenith angle. [8]
Comment on the effect of zenith angle and visibility on the samples' apparent colour. Describe the colour shifts in terms of colour names (i.e., colours tend to shift towards purple), rather than just giving coordinate numbers. At which zenith angles does the colour shift become prominent? Towards which colours do the shift converge for the different atmospheres? Explain the differences in these shifts. [4]
The picture wallpic shows a wall and some plants in sunshine and in shade. The picture was taken in mid February in Pretoria at approximately 06:45. The picture colour values are in the RGB coordinate system. The atmospheric visibility is estimated to be considerably better than 23 km, with the sun elevation angle between 10 and 20 degrees. A second picture is marked up with dots in the sunlit and shady parts of the image. Calculate the CIE colour coordinates for a number of samples in the four identified regions. Compare the colour coordinates in the four regions and comment on the shift in colour coordinates. [10]
[40]
In [29]:
# to prepare the calculation environment
import numpy as np
import pandas as pd
from scipy.interpolate import interp1d
import os
import matplotlib.pyplot as plt
import matplotlib.patches as patches
from IPython.display import display
from IPython.display import Image
from IPython.display import HTML
import pyradi.ryfiles as ryfiles
import pyradi.ryplot as ryplot
import pyradi.rymodtran as rymodtran
import pyradi.ryutils as ryutils
import pyradi.rychroma as rychroma
import pyradi.ryplanck as ryplanck
from scipy import interpolate
%matplotlib inline
%config InlineBackend.figure_format = 'png'
import PIL
#make pngs at required dpi
import matplotlib as mpl
mpl.rc("savefig", dpi=75)
mpl.rc('figure', figsize=(10,8))
doStars = False
In [30]:
# toload the sample spectral data and obtain the wavelength vector for calculations
tgzFilename = 'DP05.tgz'
destinationDir = '.'
tarFilename = 'DP05.tar'
url = 'https://raw.githubusercontent.com/NelisW/pyradi/master/pyradi/data/EOSystemAnalysisDesign-Data/'
dlNames = ryfiles.downloadUntar(tgzFilename, url, destinationDir, tarFilename)
# print('filesAvailable are {}'.format(dlNames))
wlSmpl = ryfiles.loadColumnTextFile('atmo-colour/samplesVis.txt', loadCol=[0], comment='%').reshape(-1,1)
samplesSelect = [2,6,7,8,9,10,11]
samples = ryfiles.loadColumnTextFile('atmo-colour/samplesVis.txt', loadCol=samplesSelect, comment='%')
samplesTxt=ryfiles.loadHeaderTextFile('atmo-colour/samplesVis.txt', loadCol=samplesSelect, comment='%')
# print(samplesTxt)
# print(np.hstack((wlSmpl,samples )))
# print(np.hstack((wlSmpl )))
# colourspace = 'CIE RGB'
colourspace = 'Adobe RGB'
# colourspace = 'sRGB'
Colour is an elusive property - different people perceive colour differently. Several colour standards have been proposed to minimise the effect of subjective interpretation. A commonly used standard is the CIE 1931 tristimulus colour definitions of a 'standard human observer' (shown below). In this diagram the monochromatic colours all lie on the horseshoe border.
The RGB system used in display screens are only capable of limited colour display, such as the triangle below. The colours in the following graph are not rendered correctly, because the picture was prepared for display on an RGB display that cannot display all colours. Similarly, most printers are unable to render all possible colours in the CIE 1931 system.
In [31]:
display(Image(filename="atmo-colour/CIE-planckianlocus01.png", width=400, height=400))
In the CIE 1931 colour system any colour can be expressed by two coordinates $x$ and $y$. These coordinates describe only colour, not the value (darkness/intensity). Two samples may have the same green colour, but at two different levels of brightness or darkness.
The tristimulus curves represent the colour response of the human eye. It is interesting to note that the red colour stimulus has a significant component towards the blue part of the spectrum.
The tristimulus curves are used to calculate the $x$ and $y$ colour coordinates given the spectral radiance emanating from the surface of the object.
ciebar[:,0] is the wavelength vector, ciebar[:,1] is the x-bar vector, etc.
In [32]:
# to plot the CIE tristimulus curves
cietriplt = ryplot.Plotter(1, 1,2, figsize=(12,4))
ciebarwn = rychroma.loadCIEbar(1e4/wlSmpl, stype='wn')
ciebarwl = rychroma.loadCIEbar(wlSmpl, stype='wl')
cietriplt.plot(2, ciebarwn[:,0], ciebarwn[:,1:4], "CIE tristimulus values",
r'Wavenumber cm$^{-1}$', r'Response', plotCol = ['r','g','b'],
label=['$\\bar{x}$', '$\\bar{y}$', '$\\bar{z}$'],legendAlpha=0.5,maxNX=5);
cietriplt.plot(1, ciebarwl[:,0], ciebarwl[:,1:4], "CIE tristimulus values",
r'Wavelength $\mu$m', r'Response', plotCol = ['r','g','b'],
label=['$\\bar{x}$', '$\\bar{y}$', '$\\bar{z}$'],legendAlpha=0.5);
To calculate the xyY color coordinates of a color given the spectrum, proceed as follows: \begin{eqnarray} X(T)&=&\int_0^\infty \bar{x}_\lambda E_\lambda d\lambda,\\ Y(T)&=&\int_0^\infty \bar{y}_\lambda E_\lambda d\lambda, \;\; {\rm and}\\ Z(T)&=&\int_0^\infty \bar{z}_\lambda E_\lambda d\lambda, \end{eqnarray} where $\bar{x}_\lambda$, $\bar{y}_\lambda$, and $\bar{z}_\lambda$ are the color-matching functions of the CIE standard colorimetric observer and $E$ are the spectral the Modtran irradiance values. Note that any scale factors normalises out in the division. The xyz chromaticity color coordinates can then be calculated by \begin{eqnarray} x&=&\frac{X}{X+Y+Z},\\ y&=&\frac{Y}{X+Y+Z}, \;\; {\rm and}\\ z&=&\frac{Z}{X+Y+Z} = 1-x-y, \end{eqnarray} where $x$ and $y$ define the color coordinates in the xy chromaticity diagram.
From wikipedia the conversion between RGB and CIE coordinates are given by:
The numbers in the conversion matrix below are exact, with the number of digits specified in CIE standards.[10]
\begin{equation} \begin{bmatrix}X\\Y\\Z\end{bmatrix}=\frac{1}{b_{21}} \begin{bmatrix} b_{11}&b_{12}&b_{13}\\ b_{21}&b_{22}&b_{23}\\ b_{31}&b_{32}&b_{33} \end{bmatrix} \begin{bmatrix}R\\G\\B\end{bmatrix}=\frac{1}{0.17697} \begin{bmatrix} 0.49&0.31&0.20\\ 0.17697&0.81240&0.01063\\ 0.00&0.01&0.99 \end{bmatrix} \begin{bmatrix}R\\G\\B\end{bmatrix} \end{equation}The above matrix is exactly specified in standards, but converting from RGB to XYZ is not exact:
\begin{equation} \begin{bmatrix}R\\G\\B\end{bmatrix} = \begin{bmatrix} 0.41847 & -0.15866 & -0.082835\\ -0.091169 & 0.25243 & 0.015708\\ 0.00092090 & -0.0025498 & 0.17860 \end{bmatrix} \cdot \begin{bmatrix}X\\Y\\Z\end{bmatrix} \end{equation}If the RGB coordinates of a pixel picture is known it can be converted to the CIE $(x,y)$ coordinates, using the transformation shown above. Note however, that the camera built in white balance adjustment and the spectral response of the camera detector will affect the RGB colour coordinates. Hence, this means to determine colour coordinates is inaccurate at best!
In [ ]:
In [33]:
# to plot the sRGB gamma response
rgb = np.linspace(0.,255.,500)
rgb = rgb / np.max(rgb)
# rgbo = np.where(rgb<=0.003138, 12.92*rgb, 1.055*rgb**(1./2.4)-0.055)
rgbo = np.where(rgb<=0.04045, rgb/12.92,((rgb+0.055)/1.055)**2.4)
rgb22 = rgb**2.2
p = ryplot.Plotter(1,1,1,figsize=(8,3));
p.plot(1,rgb,rgbo,'','','',plotCol='r',label=['sRGB nonlinear']);
p.plot(1,rgb,rgb22,'','','',plotCol='b',label=['$\gamma$=2.2']);
Chromaticity and value are orthogonal descriptions of a surface colour. This is best illustrated by viewing the same chromaticity at different values as in the following patch plot. Each vertical bar is has the same chromoticity but different value.
In [34]:
# to plot the samples at different values
# calculate the sample intrinsic xy colour coordinates
xii = np.zeros((samples.shape[1]))
yii = np.zeros((samples.shape[1]))
for iSmpl in range(samples.shape[1]):
# construct interpolation fn
smplfn = interp1d(1e4/wlSmpl.reshape(-1,),samples[:,iSmpl], bounds_error=False)
# get xy coords
[ xii[iSmpl], yii[iSmpl], _]=\
rychroma.chromaticityforSpectralL(ciebarwn[:,0],
smplfn(ciebarwn[:,0]).reshape(-1, 1) * ryplanck.planck(ciebarwn[:,0],6200.,'en').reshape(-1,1),
ciebarwn[:,1], ciebarwn[:,2], ciebarwn[:,3])
# convert to rgb
rgbi = rychroma.CIExy2rgb(np.hstack((xii[iSmpl].reshape(-1,1),yii[iSmpl].reshape(-1,1))),system=colourspace)
rgbi = np.where(rgbi<0,0,rgbi)
if iSmpl==0:
rgbii = rgbi
else:
rgbii = np.vstack((rgbii,rgbi))
xysamples = np.hstack((xii.reshape(-1,1),yii.reshape(-1,1) ))
rgbii /= rgbii.max(axis=0)
values = np.linspace(0.1, 1., 10 )
dsamp = 1. / samples.shape[1]
dscen = 1. / 10.
ciexyplt = ryplot.Plotter(1, 1, 1, figsize=(12,6))
pp = ciexyplt.emptyPlot(1)
for iSmpl in range(samples.shape[1]):
for k,value in enumerate(values):
rgbVal = tuple(value * rgbii[iSmpl,:])
pp.add_patch( patches.Rectangle((dsamp*iSmpl,dscen*k), dsamp, dscen,
color = rgbVal,fill=True ))
pp.text(.07+iSmpl/7., 1.01, samplesTxt[iSmpl],
horizontalalignment='center', fontsize=12)
pp.axis('off')
Out[34]:
Modtran5 was used to calculate the direct solar irradiance at sea level, and the corresponding atmospheric transmittance to the sun. The standard Modtran atmospheres selected for this investigation are:
The main differences between the above atmospheres are in are in the aerosol models used; the differences in temperature, pressure and relative humidity vertical profiles have little effect on the visual band transmittance (in the infrared bands these differences are not negligible).
The solar irradiance was calculated for each of the above scenarios at zenith angles of $\{0, 45, 60, 80, 88\}$, where zero zenith angle is vertically up and 88 degrees is near the horizon.
The data are stored nested lists. The first index into the list is the scenario index, with the second index being the zenith angle index.
tau[scenario][angle] nested list. irr[scenario][angle] nested list.ire[scenario][angle] nested list.The Modtran data are in wavenumber spectral units at 1 cm$^{-1}$ resolution --- which is too fine for these calculations. When reading in the transmittance and solar irradiance data the values are first convolved with a top-hat signal to average out the fine spectral detail from 1 cm$^{-1}$ to 40 cm$^{-1}$. The data is then converted to wavelength spectral units. Finally, the data is interpolated to the wavelength vector used in the samples and CIE tristimulus curves.
Load the Modtran data from the tape7 files: transmittance, irradiance through the atmosphere and irradiance for no atmosphere.
Note that the Modtran data are given in wavenumber terms: frequency in wavenumbers and spectral density in terms of wavenumbers. The CIE curves and sample spectral reflectance data are all in wavelength coordinates. Conversion between the two is required.
In [35]:
# to load the atmospheric data
scenarios = ['Subarctic Summer No Aerosol','Tropical Desert Aerosol','US Standard 1976 23 km Rural','Subarctic Summer 5km Urban']
angles = ['00','45','60','80','88']
tau = []
irr = []
ire = []
for i, scenario in enumerate(scenarios):
tauA = []
irrA = []
ireA = []
for j, angle in enumerate(angles):
filename = '.' + os.path.join(os.path.sep, 'atmo-colour', scenario, angle, 'tape7')
# print(filename, '{}'.format(angle))
#load the modtran tape7 file: wavenumber, transmittance, E-through-atmo, E-no-atmo
tape7 = rymodtran.loadtape7(filename, ['FREQ', 'TRANS', 'SOL_TR', 'SOLAR'] )
# convolve transmittance and irradiance from 1cm-1 to 40 cm-1, corresponds to cie increment
tape7[:,1], windowfn = ryutils.convolve(tape7[:,1], 1, 1, 40)
tape7[:,2], windowfn = ryutils.convolve(tape7[:,2], 1, 1, 40)
tape7[:,3], windowfn = ryutils.convolve(tape7[:,3], 1, 1, 40)
#convert irradiance from W/(cm2.cm-1) to W/(m2.cm-1)
tape7[:,2] = 1.0e4 * tape7[:,2]
tape7[:,3] = 1.0e4 * tape7[:,3]
#interpolate modtran data to samples spectral range
#first construct the interpolating function then call the function on modtran data
interpfunTt = interp1d(tape7[:,0], tape7[:,1], bounds_error=False, fill_value=0.0)
interpfunEt = interp1d(tape7[:,0], tape7[:,2], bounds_error=False, fill_value=0.0)
interpfunEx = interp1d(tape7[:,0], tape7[:,3], bounds_error=False, fill_value=0.0)
tauA.append(interpfunTt(ciebarwn[:,0] ))
irrA.append(interpfunEt(ciebarwn[:,0] ))
ireA.append(interpfunEx(ciebarwn[:,0] ))
tau.append(tauA)
irr.append(irrA)
ire.append(ireA)
The Modtran data represent a two-dimensional grid of scenarios and zenith angles. In the review below the data are graphically compared along each of the two dimensions to demonstrate the effect of the variation in each of the dimensions.
In [36]:
##
mT = ryplot.Plotter(1, 4, 1,"",figsize=(12,16))
for i, scenario in enumerate(scenarios):
for j, (angle, pcol) in enumerate(zip(angles,['y', 'r','b', 'g', 'k'])):
mT.plot(i+1, 1e4/ciebarwn[:,0] , tau[i][j], scenario,"Wavelength [$\mu$m]", "Transmittance",
label=['{} deg zenith'.format(angle)],plotCol=pcol, legendAlpha=0.5, pltaxis=[0.38,0.72, 0, 1])
In [37]:
##
mE = ryplot.Plotter(2, 2, 2,"",figsize=(12,10))
for i, scenario in enumerate(scenarios):
wl, Eire = ryutils.convertSpectralDensity(ciebarwn[:,0], ire[i][0],'nl')
mE.plot(i+1,wl , Eire, scenario,"Wavelength [$\mu$m]", "Irradiance W/(m$^2$.$\mu$m)",
label=['No atmosphere'],plotCol='m', legendAlpha=0.5, pltaxis=[0.38,0.72, 0, 2200])
for j, (angle, pcol) in enumerate(zip(angles,['y', 'r','b', 'g', 'k'])):
wl, Eirr = ryutils.convertSpectralDensity(ciebarwn[:,0], irr[i][j],'nl')
mE.plot(i+1,wl , Eirr, scenario,"Wavelength [$\mu$m]", "Irradiance W/(m$^2$.$\mu$m)",
label=['{} deg zenith'.format(angle)],plotCol=pcol, legendAlpha=0.5, pltaxis=[0.38,0.72, 0, 2200])
In [38]:
##
mT = ryplot.Plotter(1, 3,2,"",figsize=(12,12))
for j, angle in enumerate(angles):
for i, (scenario, pcol) in enumerate(zip(scenarios,['r', 'g','b','m'])):
mT.plot(j+1, 1e4/ciebarwn[:,0] , tau[i][j], '{} deg zenith'.format(angle), '', "Transmittance",
label=['{}'.format(scenario)],plotCol=pcol, legendAlpha=0.5, pltaxis=[0.38,0.72, 0, 1])
In [39]:
##
mE = ryplot.Plotter(2, 3, 2,"",figsize=(12,12))
for j, angle in enumerate(angles):
wl, Eire = ryutils.convertSpectralDensity(ciebarwn[:,0], ire[0][j],'nl')
mE.plot(j+1, wl , Eire, scenario,"Wavelength [$\mu$m]", "Irradiance W/(m$^2$.$\mu$m)",
label=['No atmosphere'],plotCol='m', legendAlpha=0.5, pltaxis=[0.38,0.72, 0, 2200])
for i, (scenario, pcol) in enumerate(zip(scenarios,['y', 'r','b', 'g', 'k'])):
wl, Eirr = ryutils.convertSpectralDensity(ciebarwn[:,0], irr[i][j],'nl')
mE.plot(j+1, wl , Eirr, '{} deg zenith'.format(angle), '', "Irradiance W/(m$^2$.$\mu$m)",
label=['{}'.format(scenario)],plotCol=pcol, legendAlpha=0.5, pltaxis=[0.38,0.72, 0, 2200])
It is evident that the solar zenith angle significantly affects the transmittance and solar irradiance for any given atmosphere. This effect is more pronounced at zenith angles exceeding 60 degrees.
Not surprising, the aerosol content in the atmosphere also affects the transmittance and solar irradiance even for vertical paths.
The samples illuminated by the sunlight 'filtered' by the atmosphere were measured with a spectroradiometer, illuminating the sample with a bright light at short distance. The samples are 'Lettuce', 'Skin Black', 'Skin White', 'Dark Green Leaf',u'Red Paint', 'Blue Nitrile Glove' (latex-like surgical glove), and 'White Paper'. The samples' diffuse reflection spectra The fruit samples all demonstrated considerable light propagation deeper into the fruit. The blue glove was located on a Spectralon white reference (note the considerable 'white' reflectance beyond 0.55~$\mu$m).
The wavelength vector used in the sample data is used in all subsequent spectral calculations; all other spectral data are interpolated to the values in this vector.
First calculate the samples' colour coordinates for no atmospheric attenuation and a 5900 K source temperature, i.e., the intrinsic daylight colour coordinates. These colour coordinates will serve as reference as the 'perfect' colour. In all graphs the filled circular marker represent the 'true' colour coordinate of the sample.
In the three blocks of code below,
(xi,yi) represent the intrinsic sample colour (no atmosphere, Modtran irradiance on top of the atmosphere). (Xi,Yi,Zi) represent the intrinsic sample (X,Y,Z) coordinates(no atmosphere, Modtran irradiance on top of the atmosphere). (xs,ys) represent the apparent sample colour (with atmosphere, Modtran irradiance through the atmosphere). (xm,ym) represent monochromatic colour (no atmosphere, single frequency source)
In [40]:
## to calculate the sample intrinsic xy colour coordinates
xi = np.zeros((samples.shape[1]))
yi = np.zeros((samples.shape[1]))
for iSmpl in range(samples.shape[1]):
smplfn = interp1d(1e4/wlSmpl.reshape(-1,),samples[:,iSmpl], bounds_error=False)
[ xi[iSmpl], yi[iSmpl], Y]=\
rychroma.chromaticityforSpectralL(ciebarwn[:,0],
smplfn(ciebarwn[:,0]).reshape(-1, 1) * ire[0][0].reshape(-1,1),
ciebarwn[:,1], ciebarwn[:,2], ciebarwn[:,3])
In [41]:
# to calculate the samples' intrinsic XYZ coordinates
Xi = np.zeros((samples.shape[1]))
Yi = np.zeros((samples.shape[1]))
Zi = np.zeros((samples.shape[1]))
for iSmpl in range(samples.shape[1]):
smplfn = interp1d(1e4/wlSmpl.reshape(-1,),samples[:,iSmpl], bounds_error=False)
[Xi[iSmpl], Yi[iSmpl], Zi[iSmpl]]=\
rychroma.XYZforSpectralL(ciebarwn[:,0],
smplfn(ciebarwn[:,0]).reshape(-1, 1) * ire[0][0].reshape(-1,1),
ciebarwn[:,1], ciebarwn[:,2], ciebarwn[:,3])
#perfect white
[ Xpw,YPw,ZPw] = rychroma.XYZforSpectralL(ciebarwn[:,0],
ire[0][0].reshape(-1,1),
ciebarwn[:,1], ciebarwn[:,2], ciebarwn[:,3])
Then calculate the samples' colour coordinates for the collection of scenarios and zenith angles.
In [42]:
# to calculate the samples' color coordinates through the atmosphere
xs = np.zeros((samples.shape[1],len(irr),len(irr[0])))
ys = np.zeros((samples.shape[1],len(irr),len(irr[0])))
for iSmpl in range(samples.shape[1]):
smplfn = interp1d(1e4/wlSmpl.reshape(-1,),samples[:,iSmpl], bounds_error=False)
for i, scenario in enumerate(scenarios):
for j, angle in enumerate(angles):
[ xs[iSmpl,i,j], ys[iSmpl,i,j], Y]=\
rychroma.chromaticityforSpectralL(ciebarwn[:,0],
smplfn(ciebarwn[:,0]).reshape(-1,1) * irr[i][j].reshape(-1,1),
ciebarwn[:,1], ciebarwn[:,2], ciebarwn[:,3])
In [43]:
# to plot the samples' spectral reflection
vSmpl = np.zeros((samples.shape[1]))
smpleplt = ryplot.Plotter(1, 1, 1, figsize=(8,4))
for iSmpl in range(samples.shape[1]):
# vSmpl[iSmpl] = (Xi[iSmpl] + Yi[iSmpl] + Zi[iSmpl]) / (Xpw+YPw+ZPw)
vSmpl[iSmpl] = (Yi[iSmpl]) / (YPw)
smpleplt.plot(1, wlSmpl, samples[:,iSmpl], "Sample reflectance", r'Wavelength $\mu$m',
r'Reflectance', label=['{}, value={:.2f}'.format(samplesTxt[iSmpl],vSmpl[iSmpl])],legendAlpha=0.5);
# print(samples)
Next calculate the colour coordinates for monochromatic wavelengths (the horseshoe boundary).
In [44]:
# to calculate the monochromatic line in the xy plot
xm = np.zeros(wlSmpl.shape)
ym = np.zeros(wlSmpl.shape)
#create a series of data points with unity at specific wavelength
for iWavel in range(wlSmpl.shape[0]):
monospectrum = np.zeros(wlSmpl.shape)
monospectrum[iWavel] = 1
#calc xy for single mono wavelength point
[xm[iWavel],ym[iWavel],Y] = rychroma.chromaticityforSpectralL(ciebarwn[:,0],
monospectrum, ciebarwn[:,1], ciebarwn[:,2], ciebarwn[:,3])
This function plots the colour coordinates for collections of scenarios or collections of zenith angles. Two plots are shown side by side; the left plot is the CIE xy colour coordinates for the samples and the right plot is an estimate of what the colours will look like. Please note that the estimated colour as shown here is affected by the limitations of the display or printing limitations of the RGB colour scheme. Furthermore, the xy colour coordinates do not prescribe value or scale. For the purpose of these plots, each sample was allocated an arbitrary value to demonstrate approximately the correct colour.
In [45]:
# to plot the CIE xy coordinates for all combinations of samples and atmospheres
def plotCIECoords(samples, xs, ys, scenarios, angles, scenarioSel, angleSel, vss):
ciexyplt = ryplot.Plotter(4, 1, 2, figsize=(12,6))
styleSource=['*', '<', '^', 'v','>']
styleSample=['y-','k-', 'm', 'g-','r-','b-','c-', ]
labelSample = ['A','B','C','D','E','F','G' ]
nulCC = np.asarray([0.])
zang = ['0$^\circ{}$', '45$^\circ{}$', '60$^\circ{}$', '80$^\circ{}$', '88$^\circ{}$']
for k,pltaxis in enumerate([[0, 0.9, 0, 0.9]]):
if scenarioSel >= 0:
for j, angle in enumerate(angles):
legend=['{} deg zenith'.format(angle)]
ciexyplt.plot(k+1,nulCC,nulCC,
"CIE chromaticity diagram {}".format(scenarios[scenarioSel]), r'x',r'y',
['w-'],label=legend,legendAlpha=0.5,
pltaxis=pltaxis, markers=[styleSource[j]])
if angleSel >= 0:
for i, scenario in enumerate(scenarios):
legend=['{} {} '.format(i,scenario)]
ciexyplt.plot(k+1,nulCC,nulCC,
"CIE chromaticity diagram {} deg zenith".format(angles[angleSel]), r'x',r'y',
['w-'],label=legend,legendAlpha=0.5,
pltaxis=pltaxis, markers=[styleSource[i]])
ciexyplt.plot(k+1,nulCC,nulCC,'', '','', ['w-'],
label=['True colour'],legendAlpha=0.5,
pltaxis=pltaxis, markers=['o'])
ciexyplt.resetPlotCol()
#plot monochromatic horseshoe
ciexyplt.plot(k+1, xm, ym, plotCol=['k-'],pltaxis=pltaxis)
# plot achromatic point - no colour
ciexyplt.plot(1, np.asarray([0.333, 0.333]), np.asarray([0.283, 0.383]), plotCol=['k-'])
ciexyplt.plot(1, np.asarray([0.283, 0.383]), np.asarray([0.333, 0.333]), plotCol=['k-'])
#plot chromaticity loci for samples
for iSmpl in range(samples.shape[1]):
label=samplesTxt[iSmpl]
ciexyplt.plot(k+1,xi[iSmpl],yi[iSmpl],'', '','', plotCol=[styleSample[iSmpl]],
pltaxis=pltaxis, markers=['o'])
if scenarioSel >= 0:
ciexyplt.plot(k+1,xs[iSmpl,scenarioSel,:],ys[iSmpl,scenarioSel,:],legendAlpha=0.5,
pltaxis=pltaxis,label=['{} {}'.format(labelSample[iSmpl],label)],
plotCol=[styleSample[iSmpl]])
if angleSel >= 0:
ciexyplt.plot(k+1,xs[iSmpl,:,angleSel],ys[iSmpl,:,angleSel],legendAlpha=0.5,
pltaxis=pltaxis,
label=['{} {}'.format(labelSample[iSmpl],label)], plotCol=[styleSample[iSmpl]])
#plot source markers
for iSmpl in range(samples.shape[1]):
if scenarioSel >= 0:
for j, angle in enumerate(angles):
# print(j,angle,styleSource[j])
ciexyplt.plot(k+1,xs[iSmpl,scenarioSel,j],ys[iSmpl,scenarioSel,j],
"CIE chromaticity diagram {}".format(scenarios[scenarioSel]), r'x',r'y',
[styleSample[iSmpl]],label='',
pltaxis=pltaxis, markers=[styleSource[j]])
if angleSel >= 0:
for i, scenario in enumerate(scenarios):
ciexyplt.plot(k+1,xs[iSmpl,i,angleSel],ys[iSmpl,i,angleSel],
"CIE chromaticity diagram {} deg zenith".format(angles[angleSel]), r'x',r'y',
[styleSample[iSmpl]],label='',
pltaxis=pltaxis, markers=[styleSource[i]])
pp = ciexyplt.emptyPlot(2)
for iSmpl,vs in enumerate(vss):
smplLabel=samplesTxt[iSmpl]
dsamp = 1. / samples.shape[1]
cscl = 1
rgbi = rychroma.CIExy2rgb(np.hstack((xi[iSmpl].reshape(-1,1),yi[iSmpl].reshape(-1,1))),system=colourspace)
rgbi = np.where(rgbi<0,0,rgbi)
if scenarioSel >= 0:
for i, scenario in enumerate(scenarios):
xp = xs[iSmpl,scenarioSel,:]
yp = ys[iSmpl,scenarioSel,:]
rgb = rychroma.CIExy2rgb(np.hstack((xp.reshape(-1,1),yp.reshape(-1,1))),system=colourspace)
rgb = np.where(rgb<0,0,rgb)
dangl = 1. / rgb.shape[0]
for k in range(rgb.shape[0]):
# rgbVal = tuple(cscl * vSmpl[iSmpl] * rgb[k,:])
rgbVal = tuple(cscl * vs * rgb[k,:])
pp.add_patch( patches.Rectangle((dsamp*iSmpl,dangl*k), dsamp, dangl,
color = rgbVal,fill=True ))
# zenith angles down the side of the patches
pp.text(-0.05, .1+k/5., zang[k], horizontalalignment='center', fontsize=12)
#samples across the top of the patches
pp.text(.07+iSmpl/7., 1.01, labelSample[iSmpl],
horizontalalignment='center', fontsize=12)
# pp.text(1, 1, r'11', horizontalalignment='center', fontsize=15)
# pp.text(0, 0, r'00', horizontalalignment='center', fontsize=15)
# pp.text(1, 0, r'10', horizontalalignment='center', fontsize=15)
# pp.text(0, 1, r'01', horizontalalignment='center', fontsize=15)
if angleSel >= 0:
for j, angle in enumerate(angles):
xp = xs[iSmpl,:,angleSel]
yp = ys[iSmpl,:,angleSel]
rgb = rychroma.CIExy2rgb(np.hstack((xp.reshape(-1,1),yp.reshape(-1,1))),system=colourspace)
rgb = np.where(rgb<0,0,rgb)
dscen = 1. / rgb.shape[0]
for k in range(rgb.shape[0]):
# rgbVal = tuple(cscl * vSmpl[iSmpl] * rgb[k,:])
rgbVal = tuple(cscl * vs * rgb[k,:])
pp.add_patch( patches.Rectangle((dsamp*iSmpl,dscen*k), dsamp, dscen,
color = rgbVal,fill=True ))
# zenith angles down the side of the patches
pp.text(-0.05, .1+(dscen*k), int(4.*dscen*k), horizontalalignment='center', fontsize=12)
#samples across the top of the patches
pp.text(.07+iSmpl/7., 1.01, labelSample[iSmpl],
horizontalalignment='center', fontsize=12)
pp.axis('off')
# plotCIECoords(samples=samples, xs=xs, ys=ys, scenarios=scenarios, angles=angles, scenarioSel=0,
# angleSel=-1, vss=[.8,.3,1.,.4,.7,.7,1] )
# plotCIECoords(samples=samples, xs=xs, ys=ys, scenarios=scenarios, angles=angles, scenarioSel=-1,
# angleSel=1, vss=[.8,.3,1,.4,.7,.7,1] )
In the graphs below the sample's intrinsic daylight colour (5900 K source and no atmosphere) is a circle marker to serve as a reference. The different variations of zenith angles or scenarios are indicted with star or triangle markers.
The first set of graphs plot the colour coordinates for the scenarios for the full set of zenith angles.
In [46]:
plotCIECoords(samples, xs, ys, scenarios, angles, 0, -1, vss=[.8,.3,1,.4,.7,.7,1.0] )
In [47]:
plotCIECoords(samples, xs, ys, scenarios, angles, 1, -1, vss=[.8,.3,1,.4,.7,.7,1] )
In [48]:
plotCIECoords(samples, xs, ys, scenarios, angles, 2, -1, vss=[.8,.3,1,.4,.7,.7,1] )
In [49]:
plotCIECoords(samples, xs, ys, scenarios, angles, 3, -1, vss=[.8,.3,1,.4,.7,.7,1] )
The second set of graphs plot the colour coordinates for given zenith angles for the full set of scenarios.
In [50]:
plotCIECoords(samples, xs, ys, scenarios, angles, -1, 0, vss=[.8,.3,1,.4,.7,.7,1] )
In [51]:
plotCIECoords(samples, xs, ys, scenarios, angles, -1, 1, vss=[.8,.3,1,.4,.7,.7,1] )
In [52]:
plotCIECoords(samples, xs, ys, scenarios, angles, -1, 2, vss=[.8,.3,1,.4,.7,.7,1])
In [53]:
plotCIECoords(samples, xs, ys, scenarios, angles, -1, 3, vss=[.8,.3,1,.4,.7,.7,1])
In [54]:
plotCIECoords(samples, xs, ys, scenarios, angles, -1, 4, vss=[.8,.3,1,.4,.7,.7,1])
From the colour coordinate graphs the following observations arise:
The colour coordinates shift towards yellow/red in the presence of atmospheres, even for 310 km visibility molecular scattering.
The colour coordinates shift towards red at large zenith angles (near the horizon).
The shift is small for low visibility atmospheres and at sun zenith angles below 45 deg.
The shift is progressively worse at larger zenith angles and/or lower visibility.
In the worst case of low visibility (5 km Urban) and large zenith angles the shift is so severe that all samples are essentially perceived as red in colour.
In [55]:
# to analyse and plot a single picture
def singlePictureAnalysis(picture,calpts,smplpts,colset,labels,cropwin,doAchro=True,doRGB=True):
#prepare for plotting
p = ryplot.Plotter(1,1,1,figsize=(15.7,8.5))
q = ryplot.Plotter(2,1,2,figsize=(12,6))
# plot monochromatic horseshoe
q.plot(1, xm, ym, '', '$x$', '$y$', plotCol=['k-'], pltaxis=[0.,0.9,0.,0.9])
q.plot(2, xm, ym, '',' $x$', '$y$', plotCol=['k-'], pltaxis=cropwin)
if doAchro:
# plot achromatic point - no colour
q.plot(1, np.asarray([0.333, 0.333]), np.asarray([0.283, 0.383]), plotCol=['k-'])
q.plot(1, np.asarray([0.283, 0.383]), np.asarray([0.333, 0.333]), plotCol=['k-'])
q.plot(2, np.asarray([0.333, 0.333]), np.asarray([0.283, 0.383]), plotCol=['k-'])
q.plot(2, np.asarray([0.283, 0.383]), np.asarray([0.333, 0.333]), plotCol=['k-'])
if doRGB:
#plot the RGB triangle: R, G, B
rgbTri = np.asarray([[255,0,0],[0,255,0],[0,0,255]])
xy = rychroma.rgb2CIExy(rgbTri, system=colourspace)
for i,ch in zip([0,1,2],['R','G','B']):
q.plot(1, np.asarray([xy[i,0], xy[(i+1)%3,0]]), np.asarray([xy[i,1], xy[(i+1)%3,1]]), plotCol=['k-'])
q.plot(2, np.asarray([xy[i,0], xy[(i+1)%3,0]]), np.asarray([xy[i,1], xy[(i+1)%3,1]]), plotCol=['k-'])
q.getSubPlot(1).text(xy[i,0], xy[i,1],ch)
p.showImage(1, picture, ptitle='');
# prepare for filled markers
markers = ryplot.Markers(markerfacecolor='k', marker='.')
markersxy = ryplot.Markers(markerfacecolor='k', marker='.')
if doStars:
# plot / record the markers for the calibration samples
for i,point in enumerate(calpts):
markers.add(point[0],point[1],markersize=20, markerfacecolor=colset[i],
fillstyle='full', markerfacecoloralt='k', marker='*')
xy = rychroma.rgb2CIExy(np.asarray([picture[point[1],point[0],0],
picture[point[1],point[0],1],
picture[point[1],point[0],2]]), system=colourspace)
markersxy.add(xy[0,0],xy[0,1],markersize=10, markerfacecolor=colset[i], marker='*')
# plot / record the markers for the test samples
for i,item in enumerate(smplpts):
q.plot(1,np.asarray([0.]),np.asarray([0.]),plotCol=colset[i],label=[labels[i]])
q.plot(2,np.asarray([0.3]),np.asarray([0.45]),plotCol=colset[i],label=[labels[i]])
for point in item:
markers.add(point[0],point[1],markersize=20, markerfacecolor=colset[i],
fillstyle='full', markerfacecoloralt='k')
xy = rychroma.rgb2CIExy(np.asarray([picture[point[1],point[0],0],
picture[point[1],point[0],1],
picture[point[1],point[0],2]]), system=colourspace)
markersxy.add(xy[0,0],xy[0,1],markersize=10, markerfacecolor=colset[i])
#plot the markers
markers.plot(p.getSubPlot(1))
markersxy.plot(q.getSubPlot(1))
markersxy.plot(q.getSubPlot(2))
In [56]:
fdate = '2016-06-25'
# this part of the code requires the sequence of image files and will not run without
if False:
fpath = '../../SMEOS-2016/colour-coords/data/Colour/{}/'.format(fdate)
#read excel file with solar info
xlsfilen = os.path.join(fpath,'Frame time {}.xlsx'.format(fdate))
dfSun = pd.read_excel(xlsfilen,'Solar Angle', index_col=None, na_values=['NA'],skiprows=0,header=8)
dfSun.dropna(inplace=True) # get rid of empty lines
dfSun['clock time dec'] = dfSun['clock time'].apply(lambda x: x.hour+x.minute/60.+x.second/3600.)
dfSun.sort_values(by=['clock time dec'], inplace=True)
sunElevFn = interpolate.interp1d(dfSun['clock time dec'],dfSun['altitude'])
midday = dfSun['clock time dec'].iloc[dfSun['altitude'].idxmax()]
sunrise = dfSun['clock time dec'].iloc[np.abs(dfSun['altitude'].iloc[0:360]).idxmin()]
sunset = dfSun['clock time dec'].iloc[np.abs(dfSun['altitude'].iloc[360:]).idxmin()]
print(midday)
print(sunrise)
print(sunset)
# # print(dfSun)
# print(dfSun.head())
# print(sunElevFn(12.))
In [57]:
# to load and process a sequence of images for temporal analysis
def loadimgesequence(fdate,location,notes,xlsfilename, doRGBPlot=True, doXYPlot=True, doElevPlot=False, minRGBval=2):
if not os.path.exists(xlsfilename):
# this part of the code requires the sequence of image files and will not run without
fpath = '../../SMEOS-2016/colour-coords/data/Colour/{}/'.format(fdate)
#read excel file with solar info
xlsfilen = os.path.join(fpath,'Frame time {}.xlsx'.format(fdate))
from openpyxl import load_workbook
wb = load_workbook(xlsfilen)
ws = wb.get_sheet_by_name("Solar Angle")
sunrise = ws['F2'].value
sunrise = sunrise.hour + sunrise.minute/60.
midday = ws['F4'].value
midday = midday.hour + midday.minute/60.
sunset = ws['F3'].value
sunset = sunset.hour + sunset.minute/60.
print('Sunrise {}\nSolar Noon {}\nSunset {}'.format(sunrise,midday,sunset))
dfSun = pd.read_excel(xlsfilen,'Solar Angle', index_col=None, na_values=['NA'],skiprows=0,header=8)
dfSun.dropna(inplace=True) # get rid of empty lines
dfSun['clock time dec'] = dfSun['clock time'].apply(lambda x: x.hour+x.minute/60.+x.second/3600.)
dfSun.sort_values(by=['clock time dec'], inplace=True)
sunElevFn = interpolate.interp1d(dfSun['clock time dec'],dfSun['altitude'])
# filenames = ryfiles.listFiles(fpath, patterns='06123*.png', recurse=1, return_folders=0, useRegex=False)
# filenames = ryfiles.listFiles(fpath, patterns='06394*.png', recurse=1, return_folders=0, useRegex=False)
filenames = ryfiles.listFiles(fpath, patterns='*.png', recurse=1, return_folders=0, useRegex=False)
filenames.sort()
cols = ['Zone','x','y','r','g','b','dtime','date','stime','elevation','sunrise','midday','sunset']
dfXYrgb = pd.DataFrame()
for ifi, filename in enumerate(filenames):
stime = filename.split('\\')[-1][0:6]
sdate = filename.split('\\')[-2]
dtime = (3600*float(stime[0:2]) + 60*float(stime[2:4]) + float(stime[4:6]))/3600.
if ifi%50 == 0:
print('{} -- {}'.format(dtime,filename))
picture = np.array(PIL.Image.open(filename).convert('RGB'))
for i,(item,label) in enumerate(zip(smplpts,labels)):
for point in item:
xy = rychroma.rgb2CIExy(np.asarray([picture[point[1],point[0],0],
picture[point[1],point[0],1],
picture[point[1],point[0],2]]), system=colourspace)
elev = sunElevFn(dtime)
res = [label,xy[0,0],xy[0,1],
picture[point[1],point[0],0], picture[point[1],point[0],1],
picture[point[1],point[0],2],dtime,fdate,stime,elev,
sunrise,midday,sunset]
dfXYrgb = dfXYrgb.append(pd.DataFrame([res],columns=cols),ignore_index=True)
dfXYrgb.to_excel(xlsfilename,sheet_name='rgbxy',float_format='%14.7e',na_rep=-1)
# print(dfXYrgb)
else:
dfXYrgb = pd.read_excel(xlsfilename,sheet_name='rgbxy')
dfXYrgb['x'].replace(-1.0,np.nan,inplace=True)
dfXYrgb['y'].replace(-1.0,np.nan,inplace=True)
# to process the temporal data
# remove all rows with NaN
dfXrNonan = dfXYrgb.dropna()
# print(dfXrNonan.head())
# get mean of the multiple samples
dfXrmean = pd.DataFrame()
dfXrmean['xmean'] = dfXrNonan.groupby(['Zone','dtime'])['x'].mean()
dfXrmean['ymean'] = dfXrNonan.groupby(['Zone','dtime'])['y'].mean()
dfXrmean['rmean'] = dfXrNonan.groupby(['Zone','dtime'])['r'].mean()
dfXrmean['gmean'] = dfXrNonan.groupby(['Zone','dtime'])['g'].mean()
dfXrmean['bmean'] = dfXrNonan.groupby(['Zone','dtime'])['b'].mean()
dfXrmean['elevation'] = dfXrNonan.groupby(['Zone','dtime'])['elevation'].mean()
sunrise = dfXrNonan['sunrise'].mean()
midday = dfXrNonan['midday'].mean()
sunset = dfXrNonan['sunset'].mean()
dfXrmean.reset_index(inplace=True)
print('Sunrise {}\nSolar Noon {}\nSunset {}'.format(sunrise,midday,sunset))
uZones = dfXrmean['Zone'].unique()
sZones = [u'Red', u'Green', u'Blue', u'White', u'Black']
if doRGBPlot:
# to plot the temporal RGB values
p = ryplot.Plotter(1,0.5+len(sZones)/3.,3,'Diurnal cycle {} {} {}'.format(fdate,location,notes),
figsize=(16,2.*len(sZones)))
for iz,zone in enumerate(sZones):
p.buildPlotCol(['r','g','b'])
elev = dfXrmean[dfXrmean['Zone']==zone]['elevation']
for clr,label in zip(['rmean','gmean','bmean'],['Red','Green','Blue']):
val = dfXrmean[dfXrmean['Zone']==zone][clr]
tim = dfXrmean[dfXrmean['Zone']==zone]['dtime']
p.plot(1+iz,tim,val,'{} sample'.format(zone),'Time','Value',label=[label],
pltaxis=[np.min(tim)-.5, np.max(tim)+.5,0,255],maxNX=5)
p.plot(1+iz,tim,255.*elev/90.,'{} sample'.format(zone),
label=['Sun elev, max={:.1f} deg'.format(np.max(elev))],plotCol='k',
pltaxis=[np.min(tim)-.5, np.max(tim)+.5,0,255],maxNX=5)
p.plot(1+iz,np.asarray([midday,midday]),np.asarray([0,255]),plotCol='m')
p.plot(1+iz,np.asarray([sunrise,sunrise]),np.asarray([0,255]),plotCol='c')
p.plot(1+iz,np.asarray([sunset,sunset]),np.asarray([0,255]),plotCol='c')
if doXYPlot:
# to plot the temporal xy values
q = ryplot.Plotter(2,0.5+len(sZones)/3.,3,'Diurnal cycle {} {} {}'.format(fdate,location,notes),
figsize=(16,2.*len(sZones)))
elev = dfXrmean[dfXrmean['Zone']==zone]['elevation']
for iz,zone in enumerate(sZones):
q.buildPlotCol(['m','c'])
elev = dfXrmean[dfXrmean['Zone']==zone]['elevation']
for clr,label in zip(['xmean','ymean'],['x','y']):
val = dfXrmean[dfXrmean['Zone']==zone][clr]
tim = dfXrmean[dfXrmean['Zone']==zone]['dtime']
valr = dfXrmean[dfXrmean['Zone']==zone]['rmean']
valg = dfXrmean[dfXrmean['Zone']==zone]['gmean']
valb = dfXrmean[dfXrmean['Zone']==zone]['bmean']
filt = np.logical_and(valr>minRGBval, valg>minRGBval)
filt = np.logical_and(filt>0, valb>minRGBval)
valf = val[filt]
timf = tim[filt]
q.plot(1+iz,timf,valf,'{} sample'.format(zone),'Time','Value',label=[label],
pltaxis=[np.min(tim)-.5, np.max(tim)+.5,0,1],maxNX=5)
q.plot(1+iz,tim,elev/90.,'{} sample'.format(zone),
label=['Sun elev, max={:.1f} deg'.format(np.max(elev))],plotCol='k',
pltaxis=[np.min(tim)-.5, np.max(tim)+.5,0,1],maxNX=5)
q.plot(1+iz,np.asarray([midday,midday]),np.asarray([0,1]),plotCol='g')
q.plot(1+iz,np.asarray([sunrise,sunrise]),np.asarray([0,1]),plotCol='b')
q.plot(1+iz,np.asarray([sunset,sunset]),np.asarray([0,1]),plotCol='r')
if doElevPlot:
# to plot the sun elevation angles
r = ryplot.Plotter(3,1,1,fdate,figsize=(12,3))
elev = dfXrmean[dfXrmean['Zone']=='Black']['elevation']
tim = dfXrmean[dfXrmean['Zone']=='Black']['dtime']
r.plot(1,tim,elev,'Sun elevation angle','Time','Elevation [deg]',
pltaxis=[np.min(tim)-.5, np.max(tim)+.5,-10,60],maxNX=5)
return dfXYrgb, dfXrmean
In [58]:
# to plot the temporal sequence on the xy diagram
def clktimefromdectime(dectime):
from math import floor
hours = floor(dectime)
minutes = floor((dectime*3600-hours*3600)/60)
seconds = floor(dectime*3600 - minutes*60 - hours*3600)
return('{:.0f}:{:.0f}:{:02.0f}'.format(hours, minutes,seconds))
In [59]:
#
def sequenceAnalysis(fdate,location,notes,dfXrmean,cropwin,doAchro=True,doRGB=True,reject=[0,0], minRGBval=2):
uZones = dfXrmean['Zone'].unique()
#prepare for plotting
q = ryplot.Plotter(2,1,2,'{} {} {}'.format(fdate,location,notes),figsize=(12,6))
# plot monochromatic horseshoe
q.plot(1, xm, ym, '', '$x$', '$y$', plotCol=['k-'], pltaxis=[0.,0.9,0.,0.9])
q.plot(2, xm, ym, '',' $x$', '$y$', plotCol=['k-'], pltaxis=cropwin)
if doAchro:
# plot achromatic point - no colour
q.plot(1, np.asarray([0.333, 0.333]), np.asarray([0.283, 0.383]), plotCol=['k-'])
q.plot(1, np.asarray([0.283, 0.383]), np.asarray([0.333, 0.333]), plotCol=['k-'])
q.plot(2, np.asarray([0.333, 0.333]), np.asarray([0.283, 0.383]), plotCol=['k-'])
q.plot(2, np.asarray([0.283, 0.383]), np.asarray([0.333, 0.333]), plotCol=['k-'])
if doRGB:
#plot the CIE RGB triangle: R, G, B
rgbTri = np.asarray([[255,0,0],[0,255,0],[0,0,255]])
xy = rychroma.rgb2CIExy(rgbTri, system=colourspace)
for i,ch in zip([0,1,2],['R','G','B']):
q.plot(1, np.asarray([xy[i,0], xy[(i+1)%3,0]]), np.asarray([xy[i,1], xy[(i+1)%3,1]]), plotCol=['k-'])
q.plot(2, np.asarray([xy[i,0], xy[(i+1)%3,0]]), np.asarray([xy[i,1], xy[(i+1)%3,1]]), plotCol=['k-'])
q.getSubPlot(1).text(xy[i,0], xy[i,1],ch)
plotcol = {'Red':'r','Green':'g','Blue':'b','White':'m','Black':'k'}
numdata = dfXrmean.shape[0]
dfXrmeanfilt = dfXrmean.iloc[reject[0]:numdata-reject[1]]
if 'dtime' in dfXrmean.columns:
dfXrmean.sort_values(by=['dtime'],inplace=True)
timSet = 'Time range: from {} to {}'.format(clktimefromdectime(np.min(dfXrmeanfilt['dtime'])),
clktimefromdectime(np.max(dfXrmeanfilt['dtime'])))
else:
dfXrmean.sort_values(by=['scount'],inplace=True)
timSet = 'Exposure range: from {} to {}'.format(np.min(dfXrmeanfilt['scount']),
np.max(dfXrmeanfilt['scount']))
for iz,zone in enumerate(uZones):
q.buildPlotCol(['m','c','r','g','b'])
valr = dfXrmeanfilt[dfXrmeanfilt['Zone']==zone]['rmean']
valg = dfXrmeanfilt[dfXrmeanfilt['Zone']==zone]['gmean']
valb = dfXrmeanfilt[dfXrmeanfilt['Zone']==zone]['bmean']
filt = np.logical_and(valr>minRGBval, valg>minRGBval)
filt = np.logical_and(filt>0, valb>minRGBval)
valx = dfXrmeanfilt[dfXrmeanfilt['Zone']==zone]['xmean']
valy = dfXrmeanfilt[dfXrmeanfilt['Zone']==zone]['ymean']
valxf = valx[filt]
valyf = valy[filt]
q.plot(1,valxf,valyf,timSet,label=[zone],plotCol=plotcol[zone])
q.plot(2,valxf,valyf,timSet,label=[zone],plotCol=plotcol[zone])
The analysis shown below, studies an opportunistic picture taken with a Samsung Note 4 with the standard camera software on automatic white balance correction. The picture was taken in mid February (summer) in Pretoria at approximately 06:45. The atmospheric visibility is estimated to be considerably better than 23 km, with the sun elevation angle between 10 and 20 degrees. The picture contains green foliage and a brick wall, both in the shade and in direct sunlight. The analysis shows three 'calibration' points to verify the coordinate system plus a 12 or 18 points on each of the test regions in the image.
In [60]:
# define the pixel sample coordinates
colset = ['y','g','m','r']
labels = ['Wall shadow', 'Wall sunlit','Plant shadow','Plant sunlit']
smplpts = [[[100,260],[280,260], [400,260], [540,260], [740,260], [820,260],
[100,375],[280,370], [400,370], [540,370], [740,370], [820,370],
[100,500],[280,490], [400,480], [540,480], [740,480], [820,480],],[
[1200,260],[1300,260], [1480,255], [1640,255], [1740,200], [1900,250],
[1340,470],[1210,470], [1480,360], [1640,360], [1740,360], [1900,360],
[1300,530],[1390,530], [1480,470], [1640,470], [1740,470], [1740,420],],[
[100,800],[280,800], [400,800], [540,830],
[100,950],[280,950], [400,950], [540,950],
[100,1000],[280,1000], [400,1000], [540,1000],],[
[1580,700],[1580,770],[1640,700],[1640,770],
[1580,850],[1580,925],[1640,850],[1640,925],
[1530,1000],[1530,1080],[1600,1000],[1600,1080],]]
#define the calibration coordinates to verify orientation and origin
calpts = [[10,10],[10,1200],[2260,550]]
# calib = [[100,100]]
#load the picture
wallpic = np.array(PIL.Image.open('atmo-colour/wallcolour.jpg').convert('RGB'))
print(wallpic.shape)
cropwin = [0.3,0.45,0.3,0.55]
singlePictureAnalysis(wallpic,calpts,smplpts,colset,labels,cropwin)
Even under the relatively 'good' there is considerable shift in the colour coordinates. The most striking shift is in the foliage colour from shadow to sunlight. This shift is somewhat extreme in the sense that the shaded colour is more bluish than normal and the sunlight colour shift to the right, towards yellow.
The wall colour shift is from a blue-cast yellowish pink, upwards towards yellow-green on the Planckian locus. The interesting observation here is that the colour coordinate clustering is tighter for the sunlit wall, than for the shaded wall. Visually it appears as if the blues and reds are more prominent on the shaded wall.
In [61]:
# to load and process a sequence of images for temporal analysis
def loadimgelinearitysequence(xlsfilename, lindataset,smplpts,labels, doRGBPlot=True, doXYPlot=True):
# for lindataset in ['Variable Exposure','Locked Exposure']:
if not os.path.exists(xlsfilename):
# this part of the code requires the sequence of image files and will not run without
fpath = '../../SMEOS-2016/colour-coords/data/Colour/2016-08-25/{}/'.format(lindataset)
# filenames = ryfiles.listFiles(fpath, patterns='06123*.png', recurse=1, return_folders=0, useRegex=False)
# filenames = ryfiles.listFiles(fpath, patterns='06394*.png', recurse=1, return_folders=0, useRegex=False)
filenames = ryfiles.listFiles(fpath, patterns='*.png', recurse=1, return_folders=0, useRegex=False)
print(fpath)
filenames.sort()
print(len(filenames))
cols = ['Zone','x','y','r','g','b','date','scount']
dfXYrgb = pd.DataFrame()
for ifi, filename in enumerate(filenames):
scount = filename.split('\\')[-1][6:].split('.')[0]
sdate = filename.split('\\')[-3]
if ifi%50 == 0:
print('{} -- {}'.format(scount,filename))
picture = np.array(PIL.Image.open(filename).convert('RGB'))
for i,(item,label) in enumerate(zip(smplpts,labels)):
for point in item:
xy = rychroma.rgb2CIExy(np.asarray([picture[point[1],point[0],0],
picture[point[1],point[0],1],
picture[point[1],point[0],2]]), system=colourspace)
res = [label,xy[0,0],xy[0,1],picture[point[1],point[0],0],
picture[point[1],point[0],1],
picture[point[1],point[0],2],fdate,scount]
dfXYrgb = dfXYrgb.append(pd.DataFrame([res],columns=cols),ignore_index=True)
dfXYrgb.to_excel(xlsfilename,sheet_name='rgbxy',float_format='%14.7e',na_rep=-1)
# print(dfXYrgb)
else:
dfXYrgb = pd.read_excel(xlsfilename,sheet_name='rgbxy')
dfXYrgb['x'].replace(-1.0,np.nan,inplace=True)
dfXYrgb['y'].replace(-1.0,np.nan,inplace=True)
# to process the temporal data
# remove all rows with NaN
dfXrNonan = dfXYrgb.dropna()
# get mean of the multiple samples
dfXrmean = pd.DataFrame()
dfXrmean['xmean'] = dfXrNonan.groupby(['Zone','scount'])['x'].mean()
dfXrmean['ymean'] = dfXrNonan.groupby(['Zone','scount'])['y'].mean()
dfXrmean['rmean'] = dfXrNonan.groupby(['Zone','scount'])['r'].mean()
dfXrmean['gmean'] = dfXrNonan.groupby(['Zone','scount'])['g'].mean()
dfXrmean['bmean'] = dfXrNonan.groupby(['Zone','scount'])['b'].mean()
dfXrmean.reset_index(inplace=True)
uZones = dfXrmean['Zone'].unique()
sZones = [u'Red', u'Green', u'Blue', u'White', u'Black']
if doRGBPlot:
# to plot the temporal RGB values
p = ryplot.Plotter(1,0.5+len(sZones)/3.,3,'Linearity verification {} {}'.format(fdate,lindataset),
figsize=(16,2.*len(sZones)))
for iz,zone in enumerate(sZones):
p.buildPlotCol(['r','g','b'])
for clr,label in zip(['rmean','gmean','bmean'],['Red','Green','Blue']):
val = dfXrmean[dfXrmean['Zone']==zone][clr]
tim = dfXrmean[dfXrmean['Zone']==zone]['scount']
p.plot(1+iz,tim,val,'{} sample'.format(zone),'Exposure number','Value',label=[label],
pltaxis=[np.min(tim)-.5, np.max(tim)+.5,0,255],maxNX=5)
if doXYPlot:
# to plot the temporal xy values
q = ryplot.Plotter(2,0.5+len(sZones)/3.,3,'Linearity verification {} {}'.format(fdate,lindataset),
figsize=(16,2.*len(sZones)))
for iz,zone in enumerate(sZones):
q.buildPlotCol(['m','c'])
for clr,label in zip(['xmean','ymean'],['x','y']):
val = dfXrmean[dfXrmean['Zone']==zone][clr]
tim = dfXrmean[dfXrmean['Zone']==zone]['scount']
q.plot(1+iz,tim,val,'{} sample'.format(zone),'Exposure number','Value',label=[label],
pltaxis=[np.min(tim)-.5, np.max(tim)+.5,0,1],maxNX=5)
return dfXYrgb, dfXrmean
In [62]:
# define the pixel sample coordinates
colset = ['k','y','r','g','b']
labels = ['Black','White','Red','Green','Blue']
smplpts = [
[ [700,500],[750,500],[800,500],[700,550],[750,550],[800,550],],
[ [1100,500],[1150,500],[1200,500],[1100,550],[1150,550],[1200,550],],
[ [1200,200],[1250,200],[1300,200],[1200,250],[1250,250],[1300,250],],
[ [700,200],[750,200],[800,200],[700,250],[750,250],[800,250],],
[ [1000,200],[1050,200],[1100,200],[1000,250],[1050,250],[1100,250],],
]
#define the calibration coordinates to verify orientation and origin
calpts = [[10,10],[10, 1070],[1910,10]]
#load the picture
picture = np.array(PIL.Image.open('atmo-colour/Variable-Exposure-image-0146.png').convert('RGB'))
cropwin = [0.25,0.7,0.3,0.55]
singlePictureAnalysis(picture,calpts,smplpts,colset,labels,cropwin)
In [63]:
#
dfXYrgb, dfXrmean = loadimgelinearitysequence('atmo-colour/Variable Exposure.xlsx', 'Variable Exposure', smplpts,labels, doRGBPlot=True, doXYPlot=True)
# for lindataset in ['Variable Exposure','Locked Exposure']:
In [64]:
#
cropwin = [0.2,0.75,0.15,0.6]
sequenceAnalysis('','','Linearity analysis: variable exposure',dfXrmean,cropwin,doAchro=True,doRGB=True,reject=[0,0])
In [65]:
#
cropwin = [0.2,0.75,0.15,0.6]
sequenceAnalysis('','','Linearity analysis: variable exposure',dfXrmean,cropwin,doAchro=True,doRGB=True,reject=[0,200])
In [ ]:
In [66]:
# define the pixel sample coordinates
colset = ['k','y','r','g','b']
labels = ['Black','White','Red','Green','Blue']
smplpts = [
[ [700,500],[750,500],[800,500],[700,550],[750,550],[800,550],],
[ [1100,500],[1150,500],[1200,500],[1100,550],[1150,550],[1200,550],],
[ [1250,200],[1300,200],[1350,200],[1250,250],[1300,250],[1350,250],],
[ [700,200],[750,200],[800,200],[700,250],[750,250],[800,250],],
[ [1000,200],[1050,200],[1100,200],[1000,250],[1050,250],[1100,250],],
]
#define the calibration coordinates to verify orientation and origin
calpts = [[10,10],[10, 1070],[1910,10]]
#load the picture
picture = np.array(PIL.Image.open('atmo-colour/Locked-Exposure-image-0035.png').convert('RGB'))
cropwin = [0.25,0.7,0.3,0.55]
singlePictureAnalysis(picture,calpts,smplpts,colset,labels,cropwin)
In [67]:
#
dfXYrgb, dfXrmean = loadimgelinearitysequence('atmo-colour/Locked Exposure.xlsx', 'Locked Exposure', smplpts,labels, doRGBPlot=True, doXYPlot=True)
# for lindataset in ['Variable Exposure','Locked Exposure']:
In [68]:
#
cropwin = [0.2,0.75,0.15,0.65]
sequenceAnalysis('','','Linearity analysis: locked exposure',dfXrmean,cropwin,doAchro=True,doRGB=True,reject=[0,0])
In [69]:
#
cropwin = [0.2,0.75,0.15,0.65]
sequenceAnalysis('','','Linearity analysis: locked exposure',dfXrmean,cropwin,doAchro=True,doRGB=True,reject=[0,0])
In [70]:
# define the pixel sample coordinates
colset = ['k','y','r','g','b']
labels = ['Black','White','Red','Green','Blue']
smplpts = [
[ [370,650],[420,650],[470,650],[370,700],[420,700],[470,700],],
[ [370,450],[420,450],[470,450],[370,500],[420,500],[470,500],],
[ [700,450],[750,450],[800,450],[700,500],[750,500],[800,500],],
[ [1300,450],[1350,450],[1400,450],[1300,500],[1350,500],[1400,500],],
[ [1750,450],[1800,450],[1850,450],[1750,500],[1800,500],[1850,500],],
]
#define the calibration coordinates to verify orientation and origin
calpts = [[10,10],[10, 1070],[1910,10]]
#load the picture
picture = np.array(PIL.Image.open('atmo-colour/2016-07-21-092256image-0201.png').convert('RGB'))
cropwin = [0.15,0.5,0.2,0.45]
singlePictureAnalysis(picture,calpts,smplpts,colset,labels,cropwin)
In [71]:
#
fdate = '2016-07-21'
location = 'shaded'
notes = ''
xlsfilename = 'atmo-colour/atmocolcoeff-{}-{}.xlsx'.format(fdate,location)
dfXYrgb, dfXrmean = loadimgesequence(fdate,location,notes, xlsfilename, doRGBPlot=True, doXYPlot=True)
# print(dfXYrgb.shape)
# print(dfXYrgb.iloc[dfXYrgb.shape[0]/2:dfXYrgb.shape[0]/2+10])
# print(dfXrmean.head())
# print(uZones)
In [72]:
cropwin = [0.15,0.7,0.0,0.45]
sequenceAnalysis(fdate,location,notes,dfXrmean,cropwin,doAchro=True,doRGB=True,reject=[0,0])
In [73]:
cropwin = [0.18,0.6,.15,0.425]
sequenceAnalysis(fdate,location,notes,dfXrmean,cropwin,doAchro=True,doRGB=True,reject=[80,100])
In [74]:
# define the pixel sample coordinates
colset = ['w','y','r','g','b']
labels = ['Black','White','Red','Green','Blue']
smplpts = [
[ [450,50],[500,50],[550,50],[450,100],[500,100],[550,100]],
[ [400,375],[450,375],[500,375],[400,425],[450,425],[500,425]],
[ [700,375],[750,375],[800,375],[700,425],[750,425],[800,425],],
[ [1100,350],[1150,350],[1200,350],[1100,400],[1150,400],[1200,400],],
[ [1500,350],[1550,350],[1600,350],[1500,400],[1550,400],[1600,400],],
]
#define the calibration coordinates to verify orientation and origin
calpts = [[10,10],[10, 1070],[1910,10]]
#load the picture
picture = np.array(PIL.Image.open('atmo-colour/2016-08-30-093044image-0220.png').convert('RGB'))
cropwin = [0.15,0.5,0.2,0.45]
singlePictureAnalysis(picture,calpts,smplpts,colset,labels,cropwin)
In [75]:
#
fdate = '2016-08-30'
location = 'sunlit'
notes = ''
xlsfilename = 'atmo-colour/atmocolcoeff-{}-{}.xlsx'.format(fdate,location)
dfXYrgb, dfXrmean = loadimgesequence(fdate,location,notes, xlsfilename, doRGBPlot=True, doXYPlot=True)
print(dfXYrgb.shape)
print(dfXYrgb.iloc[dfXYrgb.shape[0]/2:dfXYrgb.shape[0]/2+10])
uZones = dfXrmean['Zone'].unique()
print(dfXrmean.head())
# print(uZones)
In [76]:
cropwin = [0.15,0.7,0.1,0.45]
sequenceAnalysis(fdate,location,notes,dfXrmean,cropwin,doAchro=True,doRGB=True,reject=[0,0])
In [77]:
cropwin = [0.15,0.7,0.1,0.45]
sequenceAnalysis(fdate,location,notes,dfXrmean,cropwin,doAchro=True,doRGB=True,reject=[40,30])
In [ ]:
In [ ]:
In [ ]:
Considering only the direct sun irradiance, there is indeed a pronounced shift towards yellow/red at sun zenith angles exceeding 60 degrees. At small and intermediate zenith angles and for good visibility atmospheres the shift is towards yellow. For the combined effect of large zenith angles and low visibility atmospheres the shift is very strongly towards red.
Human experience indicate that we tend to compensate automatically for colour shift (white balance) when viewing the real world. However, when viewing a computer screen or a picture, we do not perform white balance, possibly because we do not view the object in the picture in the same context as a fully immersive real-world view.
With the assumption of visibility better than 23 km, it should be feasibly to calculate a colour correction for a given sun zenith angle. The sun zenith angle can be readily calculated from the date, time of day, and geographic location.
If the visibility is very low (i.e., 5 km) and at large zenith angles, the shift is so strong towards red that it is not feasible to correct the shift accurately. The original colours in the image are distorted so severely that manipulation would only degrade the image.
In [78]:
# you only need to do this once
#!pip install --upgrade version_information
%load_ext version_information
%version_information numpy, scipy, matplotlib, pyradi
Out[78]:
In [ ]:
In [ ]: