Add CSP modelling capabilities.

RELEASE_NOTES.rst

@@ -14,7 +14,7 @@ Release Notes
 * Atlite now supports calculating dynamic line ratings based on the IEEE-738 standard (https://github.com/PyPSA/atlite/pull/189).
 * The wind feature provided by ERA5 now also calculates the wind angle `wnd_azimuth` in range [0 - 2π) spanning the cirlce from north in clock-wise direction (0 is north, π/2 is east, -π is south, 3π/2 is west).
 * A new intersection matrix function was added, which works similarly to incidence matrix but has boolean values.
-
+* Atlite now supports two CSP (concentrated solar power) technologies, solar tower and parabolic trough. See (https://atlite.readthedocs.io/en/latest/examples/working-with-csp.html) for details.
 
 * Automated upload of code coverage reports via Codecov.
 * DataArrays returned by `.pv(...)` and `.wind(...)` now have a clearer name and 'units' attribute.

atlite/__init__.py

@@ -1,6 +1,6 @@
 # -*- coding: utf-8 -*-
 
-# SPDX-FileCopyrightText: 2016-2019 The Atlite Authors
+# SPDX-FileCopyrightText: 2016-2021 The Atlite Authors
 #
 # SPDX-License-Identifier: GPL-3.0-or-later
 
@@ -15,7 +15,7 @@
 
 from .cutout import Cutout
 from .gis import compute_indicatormatrix, regrid, ExclusionContainer
-from .resource import windturbines, solarpanels
+from .resource import windturbines, solarpanels, cspinstallations
 
 from .version import version as __version__
 

atlite/convert.py

@@ -1,6 +1,6 @@
 # -*- coding: utf-8 -*-
 
-# SPDX-FileCopyrightText: 2016-2019 The Atlite Authors
+# SPDX-FileCopyrightText: 2016-2021 The Atlite Authors
 #
 # SPDX-License-Identifier: GPL-3.0-or-later
 
@@ -33,8 +33,14 @@
 
 from . import hydro as hydrom
 from . import wind as windm
+from . import csp as cspm
 
-from .resource import get_windturbineconfig, get_solarpanelconfig, windturbine_smooth
+from .resource import (
+    get_cspinstallationconfig,
+    get_windturbineconfig,
+    get_solarpanelconfig,
+    windturbine_smooth,
+)
 
 
 def convert_and_aggregate(
@@ -502,6 +508,95 @@ def pv(cutout, panel, orientation, clearsky_model=None, **params):
     )
 
 
+# solar CSP
+def convert_csp(ds, installation):
+
+    solar_position = SolarPosition(ds)
+
+    tech = installation["technology"]
+    if tech == "parabolic trough":
+        irradiation = ds["influx_direct"]
+    elif tech == "solar tower":
+        irradiation = cspm.calculate_dni(ds, solar_position)
+    else:
+        raise ValueError(f'Unknown CSP technology option "{tech}".')
+
+    # Determine solar_position dependend efficiency for each grid cell and time step
+    efficiency = installation["efficiency"].interp(
+        altitude=solar_position["altitude"], azimuth=solar_position["azimuth"]
+    )
+
+    # Thermal system output
+    da = efficiency * irradiation
+
+    # output relative to reference irradiance
+    da /= installation["r_irradiance"]
+
+    # Limit output to max of reference irradiance
+    da = da.clip(max=1.0)
+
+    # Fill NaNs originating from DNI or solar positions outside efficiency bounds
+    da = da.fillna(0.0)
+
+    da.attrs["units"] = "kWh/kW_ref"
+    da = da.rename("specific generation")
+
+    return da
+
+
+def csp(cutout, installation, technology=None, **params):
+    """
+    Convert downward shortwave direct radiation into a csp generation time-series.
+
+    Parameters
+    ----------
+    installation: str or xr.DataArray
+        CSP installation details determining the solar field efficiency dependent on
+        the local solar position. Can be either the name of one of the standard
+        installations provided through `atlite.cspinstallationsPanel` or an
+        xarray.DataArray with 'azimuth' (in rad) and 'altitude' (in rad) coordinates
+        and an 'efficiency' (in p.u.) entry.
+    technology: str
+        Overwrite CSP technology from the installation configuration. The technology
+        affects which direct radiation is considered. Either 'parabolic trough' (DHI)
+        or 'solar tower' (DNI).
+
+    Returns
+    -------
+    csp : xr.DataArray
+        Time-series or capacity factors based on additional general
+        conversion arguments.
+
+    Note
+    ----
+    You can also specify all of the general conversion arguments
+    documented in the `convert_and_aggregate` function.
+
+    References
+    ----------
+    [1] Tobias Hirsch (ed.). SolarPACES Guideline for Bankable STE Yield Assessment,
+    IEA Technology Collaboration Programme SolarPACES, 2017.
+    URL: https://www.solarpaces.org/csp-research-tasks/task-annexes-iea/task-i-solar-thermal-electric-systems/solarpaces-guideline-for-bankable-ste-yield-assessment/
+
+    [2] Tobias Hirsch (ed.). CSPBankability Project Report, DLR, 2017.
+    URL: https://www.dlr.de/sf/en/desktopdefault.aspx/tabid-11126/19467_read-48251/
+
+    """
+
+    if isinstance(installation, (str, Path)):
+        installation = get_cspinstallationconfig(installation)
+
+    # Overwrite technology
+    if technology is not None:
+        installation["technology"] = technology
+
+    return cutout.convert_and_aggregate(
+        convert_func=convert_csp,
+        installation=installation,
+        **params,
+    )
+
+
 # hydro
 def convert_runoff(ds, weight_with_height=True):
     runoff = ds["runoff"]

atlite/csp.py

@@ -0,0 +1,59 @@
+# -*- coding: utf-8 -*-
+
+# SPDX-FileCopyrightText: 2016-2021 The Atlite Authors
+#
+# SPDX-License-Identifier: GPL-3.0-or-later
+
+"""
+Functions for use in conjunction with csp data generation.
+"""
+
+import numpy as np
+from .pv.solar_position import SolarPosition
+
+import logging
+
+logger = logging.getLogger(__name__)
+
+
+def calculate_dni(ds, solar_position=None, altitude_threshold=3.75):
+    """
+    Calculate DNI on a perpendicular plane.
+
+    Calculate the Direct Normal Irradiance (DNI) on a plane perpendicular to the solar
+    irradiance based on solar altitude and direct solar influx on a horizontal plane.
+
+    Parameters
+    ----------
+    ds : xarray.Dataset
+        Dataset containing the direct influx (influx_direct) into a horizontal plane.
+    solar_position : xarray.Dataset (optional)
+        solar_position containing a solar 'altitude' (in rad, 0 to pi/2) for the 'ds' dataset.
+        Is calculated using atlite.pv.SolarPosition if omitted.
+    altitude_threshold : float (default: 3.75 degrees)
+        Threshold for solar altitude in degrees. Values in range (0, altitude_threshold]
+        will be set to the altitude_threshold to avoid numerical issues when dividing by
+        the sine of very low solar altitude.
+        The default values '3.75 deg' corresponds to
+        the solar altitude traversed by the sun within about 15 minutes in a location with
+        maximum solar altitude of 60 deg and 10h day time.
+    """
+
+    if solar_position is None:
+        solar_position = SolarPosition(ds)
+
+    # solar altitude expected in rad, convert degrees (easier to specifcy) to match
+    altitude_threshold = np.deg2rad(altitude_threshold)
+
+    # Sanitation of altitude values:
+    # Prevent high calculated DNI values during low solar altitudes (sunset / dawn)
+    # where sin(<low altitude>) results in a very low denominator in the DNI calculation
+    altitude = solar_position["altitude"]
+    altitude = altitude.where(lambda x: x > 0, np.nan)
+    altitude = altitude.where(lambda x: x > altitude_threshold, altitude_threshold)
+
+    # Calculate DNI and remove NaNs introduced during altitude sanitation
+    # DNI is determined either by dividing by cos(azimuth) or sin(altitude)
+    dni = ds["influx_direct"] / np.sin(altitude)
+
+    return dni

atlite/cutout.py

@@ -46,6 +46,7 @@
     temperature,
     wind,
     pv,
+    csp,
     runoff,
     solar_thermal,
     soil_temperature,
@@ -629,6 +630,8 @@ def layout_from_capacity_list(self, data, col="Capacity"):
 
     pv = pv
 
+    csp = csp
+
     runoff = runoff
 
     hydro = hydro

atlite/pv/solar_position.py

@@ -1,6 +1,6 @@
 # -*- coding: utf-8 -*-
 
-# SPDX-FileCopyrightText: 2016-2019 The Atlite Authors
+# SPDX-FileCopyrightText: 2016-2021 The Atlite Authors
 #
 # SPDX-License-Identifier: GPL-3.0-or-later
 
@@ -91,4 +91,5 @@ def SolarPosition(ds):
 
     vars = {da.name: da for da in [alt, az, atmospheric_insolation]}
     solar_position = xr.Dataset(vars)
+
     return solar_position

atlite/resource.py

@@ -1,6 +1,6 @@
 # -*- coding: utf-8 -*-
 
-# SPDX-FileCopyrightText: 2016-2019 The Atlite Authors
+# SPDX-FileCopyrightText: 2016-2021 The Atlite Authors
 #
 # SPDX-License-Identifier: GPL-3.0-or-later
 
@@ -28,6 +28,7 @@
 RESOURCE_DIRECTORY = Path(pkg_resources.resource_filename(__name__, "resources"))
 WINDTURBINE_DIRECTORY = RESOURCE_DIRECTORY / "windturbine"
 SOLARPANEL_DIRECTORY = RESOURCE_DIRECTORY / "solarpanel"
+CSPINSTALLATION_DIRECTORY = RESOURCE_DIRECTORY / "cspinstallation"
 
 
 def get_windturbineconfig(turbine):
@@ -83,6 +84,61 @@ def get_solarpanelconfig(panel):
     return conf
 
 
+def get_cspinstallationconfig(installation):
+    """Load the 'installation'.yaml file from local disk to provide the system efficiencies.
+
+    Parameters
+    ----------
+    installation : str
+        Name of CSP installation kind. Must correspond to name of one of the files
+        in resources/cspinstallation.
+
+    Returns
+    -------
+    config : dict
+        Config with details on the CSP installation.
+    """
+
+    if isinstance(installation, str):
+        if not installation.endswith(".yaml"):
+            installation += ".yaml"
+
+    installation = CSPINSTALLATION_DIRECTORY / installation
+
+    # Load and set expected index columns
+    with open(installation, "r") as f:
+        config = yaml.safe_load(f)
+
+    config["path"] = installation
+
+    ## Convert efficiency dict to xr.DataArray and convert units to deg -> rad, % -> p.u.
+    da = pd.DataFrame(config["efficiency"]).set_index(["altitude", "azimuth"])
+
+    # Handle as xarray DataArray early - da will be 'return'-ed
+    da = da.to_xarray()["value"]
+
+    # Solar altitude + azimuth expected in deg for better readibility
+    # calculations use solar position in rad
+    # Convert da to new coordinates and drop old
+    da = da.rename({"azimuth": "azimuth [deg]", "altitude": "altitude [deg]"})
+    da = da.assign_coords(
+        {
+            "altitude": np.deg2rad(da["altitude [deg]"]),
+            "azimuth": np.deg2rad(da["azimuth [deg]"]),
+        }
+    )
+    da = da.swap_dims({"altitude [deg]": "altitude", "azimuth [deg]": "azimuth"})
+
+    da = da.chunk("auto")
+
+    # Efficiency unit from % to p.u.
+    da /= 1.0e2
+
+    config["efficiency"] = da
+
+    return config
+
+
 def solarpanel_rated_capacity_per_unit(panel):
     # unit is m^2 here
 
@@ -319,3 +375,6 @@ def get_oedb_windturbineconfig(search=None, **search_params):
 _oedb_turbines = None
 windturbines = arrowdict({p.stem: p for p in WINDTURBINE_DIRECTORY.glob("*.yaml")})
 solarpanels = arrowdict({p.stem: p for p in SOLARPANEL_DIRECTORY.glob("*.yaml")})
+cspinstallations = arrowdict(
+    {p.stem: p for p in CSPINSTALLATION_DIRECTORY.glob("*.yaml")}
+)