start developing sky_map.py

move sky_grid and show_sky in sky_map

src/alinea/astk/meteorology/sky_irradiance.py

@@ -316,11 +316,11 @@ def actual_sky_irradiances(dates=None, daydate=_daydate, ghi=None,
     return df.loc[:, ('ghi', 'dhi', 'dni')]
 
 
-def sky_irradiances(dates=None, daydate=_daydate, ghi=None, dhi=None, ppfd=None,
-                           attenuation=None,
-                           pressure=101325, temp_dew=None, longitude=_longitude,
-                           latitude=_latitude, altitude=_altitude,
-                           timezone=_timezone, with_pvlib=True):
+def sky_irradiance(dates=None, daydate=_daydate, ghi=None, dhi=None, ppfd=None,
+                   attenuation=None,
+                   pressure=101325, temp_dew=None, longitude=_longitude,
+                   latitude=_latitude, altitude=_altitude,
+                   timezone=_timezone, with_pvlib=True):
     """ Estimate variables related to sky irradiance.
 
     Args:

src/alinea/astk/meteorology/sky_luminance.py

@@ -13,12 +13,9 @@
 # ==============================================================================
 """ A collection of equation for modelling distribution of sky luminance
 """
-from __future__ import division
-
 import numpy
-from matplotlib import pyplot as plt
-
-from alinea.astk.meteorology.sky_irradiance import horizontal_irradiance, all_weather_sky_clearness, f_clear_sky, all_weather_sky_brightness
+from alinea.astk.meteorology.sky_irradiance import (horizontal_irradiance, all_weather_sky_clearness, f_clear_sky,
+                                                    all_weather_sky_brightness)
 
 
 def cie_luminance_gradation(z, a=4, b=-0.7):
@@ -37,28 +34,6 @@ def _f(z):
     return _f(z) / _f(0)
 
 
-def _ksi(sky_zenith, sky_azimuth, sun_zenith=0, sun_azimuth=0):
-    """acute angle between an array of sky elements and the sun """
-
-    def _xyz(zenith, azimuth):
-        theta = numpy.radians(zenith)
-        phi = numpy.radians(azimuth)
-        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, axis=2) * numpy.linalg.norm(v2))
-        angle = numpy.arccos(numpy.clip(x, -1, 1))
-        return numpy.degrees(angle)
-
-    v_sky = numpy.stack(_xyz(sky_zenith, sky_azimuth), axis=2)
-    v_sun = _xyz(sun_zenith, sun_azimuth)
-
-    return _angle3(v_sky, v_sun)
-
-
 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
@@ -78,7 +53,7 @@ def _f(k):
     return _f(ksi) / _f(ksi_sun)
 
 
-def cie_relative_luminance(sky_zenith, sky_azimuth=None, sun_zenith=None,
+def cie_relative_luminance(sky_zenith=None, grid=None, sun_zenith=None,
                            sun_azimuth=None, type='soc'):
     """ cie relative luminance of a sky element relative to the luminance
     at zenith
@@ -91,9 +66,13 @@ def cie_relative_luminance(sky_zenith, sky_azimuth=None, sun_zenith=None,
     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 sky_zenith is None and grid is None:
+        raise ValueError('Either sky_zenith or grid should be passed')
+    sky_azimuth = None
+    if grid is not None:
+        _, _, _, sky_zenith = grid
+    else:
+        sky_zenith = numpy.array(sky_zenith)
 
     if type == 'soc':
         ab = {'a': 4, 'b': -0.7}
@@ -104,13 +83,13 @@ def cie_relative_luminance(sky_zenith, sky_azimuth=None, sun_zenith=None,
     else:
         raise ValueError('unknown sky_type:' + type)
 
-    gradation = cie_luminance_gradation(numpy.array(sky_zenith), **ab)
+    gradation = cie_luminance_gradation(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)
+        ksi = ksi_grid(grid, sun_zenith, sun_azimuth)
         indicatrix = cie_scattering_indicatrix(ksi, ksi_sun=ksi_sun, **cde)
 
     return gradation * indicatrix
@@ -180,13 +159,35 @@ def _awfit(p1, p2, p3, p4, zen, br):
     return a, b, c, d, e
 
 
-def all_weather_relative_luminance(sky_zenith, sky_azimuth, sun_zenith, sun_azimuth, clearness, brightness):
+def ksi_grid(grid, sun_zenith=0, sun_azimuth=0):
+    """acute angle between an array of sky elements and the sun """
+
+    def _xyz(zenith, azimuth):
+        theta = numpy.radians(zenith)
+        phi = numpy.radians(azimuth)
+        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, axis=2) * numpy.linalg.norm(v2))
+        angle = numpy.arccos(numpy.clip(x, -1, 1))
+        return numpy.degrees(angle)
+
+    _, _, sky_azimuth, sky_zenith = grid
+    v_sky = numpy.stack(_xyz(sky_zenith, sky_azimuth), axis=2)
+    v_sun = _xyz(sun_zenith, sun_azimuth)
+
+    return _angle3(v_sky, v_sun)
+
+
+def all_weather_relative_luminance(grid, sun_zenith, sun_azimuth, clearness, brightness):
     """All weather relative luminance of a sky element relative to the luminance
     at zenith
 
     Args:
-        sky_zenith : zenith angle of the sky element (deg)
-        sky_azimuth: azimuth angle of the sky element (deg)
+        grid: a (azimuth, zenith, az_c, z_c) tuple, such as returned by astk.sky_grid
         sun_zenith : zenith angle of the sun (deg)
         sun_azimuth: azimuth angle of the sun (deg)
         clearness: sky clearness as defined in Perez et al. (1993
@@ -197,38 +198,26 @@ def all_weather_relative_luminance(sky_zenith, sky_azimuth, sun_zenith, sun_azim
             validation", Solar Energy, Volume 50, Issue 3, 1993, Pages 235-245,
 
     """
-
     a, b, c, d, e = all_weather_abcde(sun_zenith, clearness, brightness)
-    gradation = cie_luminance_gradation(numpy.array(sky_zenith), a=a, b=b)
+    _, _, _, sky_zenith = grid
+    gradation = cie_luminance_gradation(sky_zenith, a=a, b=b)
     ksi_sun = sun_zenith
-    ksi = _ksi(sky_zenith, sky_azimuth, sun_zenith, sun_azimuth)
+    ksi = ksi_grid(grid, sun_zenith, sun_azimuth)
     indicatrix = cie_scattering_indicatrix(ksi, ksi_sun=ksi_sun, c=c, d=d, e=e)
 
     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, da=1):
+def sky_luminance(grid, 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', 'all_weather')
         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'
+        astk.sky_irradiance.sky_irradiance. Needed for all sky_types except 'uoc' and 'soc'
     """
 
-    azimuth, zenith, az_c, z_c = sky_grid(daz=da, dz=da)
+    azimuth, zenith, az_c, z_c = grid
     lum = numpy.zeros_like(az_c)
     hi = 0
     irrad = sky_irradiance.copy()
@@ -239,7 +228,7 @@ def sky_luminance(sky_type='soc', sky_irradiance=None, da=1):
         t = 'soc'
         if sky_type == 'uoc':
             t = 'uoc'
-        lum = cie_relative_luminance(z_c, az_c, type=t)
+        lum = cie_relative_luminance(grid=grid, type=t)
         lum /= lum.sum()
         hi = horizontal_irradiance(lum, 90 - z_c).sum()
 
@@ -250,7 +239,8 @@ def sky_luminance(sky_type='soc', sky_irradiance=None, da=1):
             # scale diffuse part to sum of diffuse irradiance
             lum *= (irrad.dhi.sum() / hi)
             for row in irrad.itertuples():
-                i, j = int(row.sun_zenith // da), int(row.sun_azimuth // da)
+                daz, dz = map(lambda x: numpy.diff(x)[0], (azimuth, zenith))
+                i, j = int(row.sun_zenith // dz), int(row.sun_azimuth // daz)
                 lum[i, j] += row.dni
             lum /= lum.sum()
             hi = horizontal_irradiance(lum, 90 - z_c).sum()
@@ -264,7 +254,7 @@ def sky_luminance(sky_type='soc', sky_irradiance=None, da=1):
                 if sky_type == 'all_weather':
                     brightness = all_weather_sky_brightness(row.dates, row.dhi, row.sun_zenith)
                     clearness = all_weather_sky_clearness(row.dni, row.dhi, row.sun_zenith)
-                    _lum = all_weather_relative_luminance(z_c, az_c,
+                    _lum = all_weather_relative_luminance(grid,
                                                           sun_zenith=row.sun_zenith,
                                                           sun_azimuth=row.sun_azimuth,
                                                           brightness=brightness,
@@ -273,7 +263,7 @@ def sky_luminance(sky_type='soc', sky_irradiance=None, da=1):
                     _hi = horizontal_irradiance(_lum, 90 - z_c).sum()
                     lum += (row.ghi / _hi * _lum)
                 else:
-                    _lum = cie_relative_luminance(z_c, az_c,
+                    _lum = cie_relative_luminance(grid=grid,
                                                    sun_zenith=row.sun_zenith,
                                                    sun_azimuth=row.sun_azimuth,
                                                    type='clear_sky')
@@ -290,21 +280,7 @@ def sky_luminance(sky_type='soc', sky_irradiance=None, da=1):
             lum /= lum.sum()
             hi = horizontal_irradiance(lum, 90 - z_c).sum()
 
-    return azimuth, zenith, lum, hi
-
-
-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 = z
-    fig, ax = plt.subplots(subplot_kw={'projection': 'polar'})
-    ax.pcolormesh(theta, r, lum, edgecolors='face', cmap=cmap, shading=shading)
-    plt.show()
+    return lum, hi
 
 
 

src/alinea/astk/sky_map.py

@@ -0,0 +1,53 @@
+# -*- 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>
+#
+# ==============================================================================
+"""Creation and plotting of sky maps
+"""
+import numpy
+from matplotlib import pyplot as plt
+
+
+def sky_grid(d_az=1, d_z=1, n_az=None, n_z=None):
+    """Azimuth and zenital grid coordinates"""
+    def _c(x):
+        return (x[:-1] + x[1:]) / 2
+    # grid cell boundaries
+    if n_az is None:
+        azimuth = numpy.linspace(0, 360, 360 // d_az + 1)
+    else:
+        azimuth = numpy.linspace(0, 360, n_az + 1)
+    if n_z is None:
+        zenith = numpy.linspace(0, 90, 90 // d_z + 1)
+    else:
+        zenith = numpy.linspace(0, 90, n_z + 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 show_sky(grid, sky, cmap='jet', shading='flat'):
+    """Display sky luminance polar image
+
+    Args:
+        grid: a (azimuth, zenith, az_c, z_c) tuple, such as returned by sky_grid
+        sky: a 2-D array of values to be plotted on the grid
+    """
+    az, z, _, _ = grid
+    theta = numpy.radians(az)
+    r = z
+    fig, ax = plt.subplots(subplot_kw={'projection': 'polar'})
+    ax.pcolormesh(theta, r, sky, edgecolors='face', cmap=cmap, shading=shading)
+    plt.show()
+
+

src/alinea/astk/sun_and_sky.py

@@ -19,9 +19,9 @@
 import numpy
 import pandas
 from alinea.astk.meteorology.sky_irradiance import (
-    sky_irradiances,
+    sky_irradiance,
     clear_sky_irradiances,
-    horizontal_irradiance, f_clear_sky)
+    horizontal_irradiance)
 from alinea.astk.meteorology.sky_luminance import cie_relative_luminance
 from alinea.astk.meteorology.sun_position import sun_position
 
@@ -297,11 +297,11 @@ def sun_sky_sources(ghi=None, dhi=None, attenuation=None, model='blended',
         Leicester, UK, 2000.
     """
 
-    sky_irr = sky_irradiances(dates=dates, daydate=daydate, ghi=ghi, dhi=dhi,
-                              attenuation=attenuation, pressure=pressure,
-                              temp_dew=temp_dew, longitude=longitude,
-                              latitude=latitude, altitude=altitude,
-                              timezone=timezone)
+    sky_irr = sky_irradiance(dates=dates, daydate=daydate, ghi=ghi, dhi=dhi,
+                             attenuation=attenuation, pressure=pressure,
+                             temp_dew=temp_dew, longitude=longitude,
+                             latitude=latitude, altitude=altitude,
+                             timezone=timezone)
 
     f_sun = sun_fraction(sky_irr)
     sun = sun_sources(irradiance=f_sun, dates=dates,

test/test_sky_irradiance.py

@@ -1,8 +1,8 @@
 import numpy
 from alinea.astk.meteorology.sky_irradiance import (
-    clear_sky_irradiances, 
-    actual_sky_irradiances, 
-    sky_irradiances)
+    clear_sky_irradiances,
+    actual_sky_irradiances,
+    sky_irradiance)
 
 
 
@@ -30,11 +30,11 @@ def test_actual_sky_irradiance():
 
 
 def test_sky_irradiances():
-    df = sky_irradiances()
+    df = sky_irradiance()
     assert len(df) == 15
     assert numpy.isclose((df.ghi/df.ppfd).mean(), 0.47, atol=0.01)
     assert df.dhi.sum() / df.ghi.sum() < 0.25
-    df = sky_irradiances(attenuation=0.2)
+    df = sky_irradiance(attenuation=0.2)
     assert df.dhi.sum() / df.ghi.sum() > 0.99
-    df2 = sky_irradiances(with_pvlib=False)
+    df2 = sky_irradiance(with_pvlib=False)
     assert len(df2) == 15