create sky_luminance for cie skies with polar plot of resulting lumin…

…ances

src/alinea/astk/Weather.py

@@ -9,7 +9,7 @@
 from __future__ import print_function
 import pandas
 import pytz
-from datetime import datetime, timedelta
+from datetime import timedelta
 
 
 from alinea.astk.TimeControl import *

src/alinea/astk/meteorology/sky_irradiance.py

@@ -47,14 +47,14 @@
 
 
 def horizontal_irradiance(normal_irradiance, elevation):
-    """ irradiance measured on an horizontal surface from a source
+    """ irradiance measured on a horizontal surface from a source
     with known elevation (degrees) and known normal irradiance
     """
     return normal_irradiance * numpy.sin(numpy.radians(elevation))
 
 
 def normal_irradiance(horizontal_irradiance, elevation):
-    """ irradiance measured on an surface perpendicular
+    """ irradiance measured on a surface perpendicular
     to a source with known elevation (degrees) and horizontal irradiance
     """
     return horizontal_irradiance / numpy.sin(numpy.radians(elevation))

src/alinea/astk/meteorology/sky_luminance.py

@@ -0,0 +1,166 @@
+# -*- python -*-
+#
+#       Copyright 2016 INRIA - CIRAD - INRA
+#
+#       Distributed under the Cecill-C License.
+#       See accompanying file LICENSE.txt or copy at
+#           http://www.cecill.info/licences/Licence_CeCILL-C_V1-en.html
+#
+#       WebSite : https://github.com/openalea-incubator/astk
+#
+#       File author(s): Christian Fournier <Christian.Fournier@supagro.inra.fr>
+#
+# ==============================================================================
+""" A collection of equation for modelling distribution of sky luminance
+"""
+from __future__ import division
+
+import numpy
+from matplotlib import pyplot as plt
+from functools import reduce
+
+from alinea.astk.meteorology.sky_irradiance import horizontal_irradiance
+
+def cie_luminance_gradation(z, a=4, b=-0.7):
+    """ function giving the dependence of the luminance of a sky element
+    to its zenith angle
+
+    CIE, 2002, Spatial distribution of daylight CIE standard general sky,
+    CIE standard, CIE Central Bureau, Vienna
+
+    z : zenith angle of the sky element (deg)
+    a, b : coefficient for the type of sky
+    """
+    def _f(z):
+        return 1 + numpy.where(z == 90, 0,
+                               a * numpy.exp(b / numpy.cos(numpy.radians(z))))
+    return _f(z) / _f(0)
+
+
+def _ksi(sky_zenith, sky_azimuth, sun_zenith=0, sun_azimuth=0):
+    """acute angle between a sky element and the sun """
+
+    def _xyz(theta, phi):
+        return (numpy.sin(theta) * numpy.cos(phi),
+                numpy.sin(theta) * numpy.sin(phi),
+                numpy.cos(theta))
+
+    def _angle3(v1, v2):
+        """acute angle between 2 3d vectors"""
+        x = numpy.dot(v1, v2) / (numpy.linalg.norm(v1) * numpy.linalg.norm(v2))
+        angle = numpy.arccos(numpy.clip(x, -1, 1))
+        return numpy.degrees(angle)
+
+    z_sun, az_sun, z_sky, az_sky = map(numpy.radians, (sun_zenith, sun_azimuth, sky_zenith, sky_azimuth))
+    v_sun = _xyz(z_sun, az_sun)
+    v_sky = _xyz(z_sky, az_sky)
+    return numpy.apply_along_axis(lambda x: _angle3(x,v_sun),0,numpy.stack(v_sky))
+
+
+def cie_scattering_indicatrix(ksi, ksi_sun=0, c=10, d=-3, e=-0.45):
+    """ function giving the dependence of the luminance
+    to its angular distance to the sun
+
+    CIE, 2002, Spatial distribution of daylight CIE standard general sky,
+    CIE standard, CIE Central Bureau, Vienna
+
+    ksi: angular distance to the sun (deg)
+    ksi_sun: angular distance of zenith to the sun, ie zenith angle of the sun (deg)
+    c, d, e : coefficient for the type of sky
+    """
+    def _f(k):
+        ksi_r = numpy.radians(k)
+        return 1 + c * (
+            numpy.exp(d * ksi_r) - numpy.exp(d * numpy.pi / 2)) + e * numpy.power(
+        numpy.cos(ksi_r), 2)
+    return _f(ksi) / _f(ksi_sun)
+
+
+def cie_relative_luminance(sky_zenith, sky_azimuth=None, sun_zenith=None,
+                           sun_azimuth=None, type='soc'):
+    """ cie relative luminance of a sky element relative to the luminance
+    at zenith
+
+    sky_zenith : zenith angle of the sky element (deg)
+    sky_azimuth: azimuth angle of the sky element (deg)
+    sun_zenith : zenith angle of the sun (deg)
+    sun_azimuth: azimuth angle of the sun (deg)
+    type is one of 'soc' (standard overcast sky), 'uoc' (uniform radiance)
+    or 'clear_sky' (standard clear sky low turbidity)
+    """
+
+    if type == 'clear_sky' and (
+                sun_zenith is None or sun_azimuth is None or sky_azimuth is None):
+        raise ValueError('Clear sky requires sun position')
+
+    if type == 'soc':
+        ab = {'a': 4, 'b': -0.7}
+    elif type == 'uoc':
+        ab = {'a': 0, 'b': -1}
+    elif type == 'clear_sky':
+        ab = {'a': -1, 'b': -0.32}
+    else:
+        raise ValueError('unknown sky_type:' + type)
+
+    gradation = cie_luminance_gradation(numpy.array(sky_zenith), **ab)
+
+    indicatrix = 1
+    if type == 'clear_sky':
+        cde = {'c': 10, 'd': -3, 'e': 0.45}
+        ksi_sun = sun_zenith
+        ksi = _ksi(sky_zenith, sky_azimuth, sun_zenith, sun_azimuth)
+        indicatrix = cie_scattering_indicatrix(ksi, ksi_sun=ksi_sun, **cde)
+
+    return gradation * indicatrix
+
+
+def sky_grid(daz=1,dz=1):
+    """Azimuth and zenital grid coordinates"""
+    def _c(x):
+        return (x[:-1] + x[1:]) / 2
+    # grid cell boundaries
+    azimuth = numpy.linspace(0, 360, 360 // daz + 1)
+    zenith = numpy.linspace(0, 90, 90 // dz + 1)
+    # grid cell centers positioned in a 2D grid matrix
+    az_c, z_c = numpy.meshgrid(_c(azimuth), _c(zenith))
+    return azimuth, zenith , az_c, z_c
+
+
+def sky_luminance(sky_type='soc', sky_irradiance=None):
+    """Sky luminance map for different conditions and sky types
+
+    Args:
+        sky_type (str): sky type, one of ('soc', 'uoc', 'clear_sky', 'sun_soc', 'blended')
+        sky_irradiance: a datetime indexed dataframe specifying sky irradiances for the period, such as returned by
+        astk.sky_irradiance.sky_irradiances. Needed for all sky_types except 'uoc' and 'soc'
+    """
+
+    azimuth, zenith, az_c, z_c = sky_grid()
+    lum = None
+    if sky_type in ('soc', 'uoc'):
+        lum = cie_relative_luminance(z_c, az_c, type=sky_type)
+    elif sky_type == 'clear_sky':
+        assert sky_irradiance is not None
+        lums = []
+        for row in sky_irradiance.itertuples():
+            lums.append(cie_relative_luminance(z_c, az_c,
+                                               sun_zenith=row.sun_zenith,
+                                               sun_azimuth=row.sun_azimuth,
+                                               type=sky_type))
+        lum = reduce(lambda a, b: a + b, lums)
+    lum /= lum.sum()
+    return azimuth, zenith, lum
+
+
+def show_sky(sky, cmap='jet', shading='flat'):
+    """Display sky luminance polar image
+
+    Args:
+        sky: a azimuth, zenith, luminance tuple
+    """
+    az, z, lum = sky
+    theta = numpy.radians(az)
+    r = numpy.sin(numpy.radians(z)) * 90
+    fig, ax = plt.subplots(subplot_kw={'projection': 'polar'})
+    ax.pcolormesh(theta, r, lum, edgecolors='face', cmap=cmap, shading=shading)
+    plt.show()

src/alinea/astk/sun_and_sky.py

@@ -23,6 +23,7 @@
     clear_sky_irradiances,
     horizontal_irradiance,
     all_weather_sky_clearness)
+from alinea.astk.meteorology.sky_luminance import cie_relative_luminance
 from alinea.astk.meteorology.sun_position import sun_position
 from six.moves import map
 
@@ -35,77 +36,6 @@
 
 
 # sky models / equations
-def cie_luminance_gradation(sky_elevation, a, b):
-    """ function giving the dependence of the luminance of a sky element
-    to its elevation angle
-
-    CIE, 2002, Spatial distribution of daylight CIE standard general sky,
-    CIE standard, CIE Central Bureau, Vienna
-
-    elevation : elevation angle of the sky element (rad)
-    a, b : coefficient for the type of sky
-    """
-    z = numpy.pi / 2 - numpy.array(sky_elevation)
-    phi_0 = 1 + a * numpy.exp(b)
-    phi_z = numpy.where(sky_elevation == 0, 1,
-                        1 + a * numpy.exp(b / numpy.cos(z)))
-    return phi_z / phi_0
-
-
-def cie_scattering_indicatrix(sun_azimuth, sun_elevation, sky_azimuth,
-                              sky_elevation, c, d, e):
-    """ function giving the dependence of the luminance
-    to its azimuth distance to the sun
-
-    CIE, 2002, Spatial distribution of daylight CIE standard general sky,
-    CIE standard, CIE Central Bureau, Vienna
-
-    elevation : elevation angle of the sky element (rad)
-    d, e : coefficient for the type of sky
-    """
-    z = numpy.pi / 2 - numpy.array(sky_elevation)
-    zs = numpy.pi / 2 - numpy.array(sun_elevation)
-    alpha = numpy.array(sky_azimuth)
-    alpha_s = numpy.array(sun_azimuth)
-    ksi = numpy.arccos(
-        numpy.cos(zs) * numpy.cos(z) + numpy.sin(zs) * numpy.sin(z) * numpy.cos(
-            numpy.abs(alpha - alpha_s)))
-
-    f_ksi = 1 + c * (
-    numpy.exp(d * ksi) - numpy.exp(d * numpy.pi / 2)) + e * numpy.power(
-        numpy.cos(ksi), 2)
-    f_zs = 1 + c * (
-    numpy.exp(d * zs) - numpy.exp(d * numpy.pi / 2)) + e * numpy.power(
-        numpy.cos(zs), 2)
-
-    return f_ksi / f_zs
-
-
-def cie_relative_luminance(sky_elevation, sky_azimuth=None, sun_elevation=None,
-                           sun_azimuth=None, type='soc'):
-    """ cie relative luminance of a sky element relative to the luminance
-    at zenith
-
-    angle in radians
-    type is one of 'soc' (standard overcast sky), 'uoc' (uniform radiance)
-    or 'clear_sky' (standard clear sky low turbidity)
-    """
-
-    if type == 'clear_sky' and (
-                sun_elevation is None or sun_azimuth is None or sky_azimuth is None):
-        raise ValueError('Clear sky requires sun position')
-
-    if type == 'soc':
-        return cie_luminance_gradation(sky_elevation, 4, -0.7)
-    elif type == 'uoc':
-        return cie_luminance_gradation(sky_elevation, 0, -1)
-    elif type == 'clear_sky':
-        return cie_luminance_gradation(sky_elevation, -1,
-                                       -0.32) * cie_scattering_indicatrix(
-            sun_azimuth, sun_elevation, sky_azimuth, sky_elevation, 10, -3,
-            0.45)
-    else:
-        raise ValueError('Unknown sky type')
 
 
 def sky_discretisation(turtle_sectors=46, nb_az=None, nb_el=None):