Dynamic line rating

RELEASE_NOTES.rst

@@ -11,6 +11,11 @@ Release Notes
 .. Upcoming Release
 .. =================
 
+* 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.
+
+
 * Automated upload of code coverage reports via Codecov.
 
 Version 0.2.5 

atlite/convert.py

@@ -15,6 +15,7 @@
 import datetime as dt
 from operator import itemgetter
 from pathlib import Path
+from dask import delayed, compute
 from dask.diagnostics import ProgressBar
 from scipy.sparse import csr_matrix
 
@@ -615,3 +616,208 @@ def hydro(
     return hydrom.shift_and_aggregate_runoff_for_plants(
         basins, runoff, flowspeed, show_progress
     )
+
+
+def convert_line_rating(
+    ds, psi, R, D=0.028, Ts=373, epsilon=0.6, alpha=0.6, per_unit=False
+):
+    """
+    Convert the cutout data to dynamic line rating time series.
+
+    The formulation is based on:
+
+    [1]“IEEE Std 738™-2012 (Revision of IEEE Std 738-2006/Incorporates IEEE Std
+        738-2012/Cor 1-2013), IEEE Standard for Calculating the Current-Temperature
+        Relationship of Bare Overhead Conductors,” p. 72.
+
+    The following simplifications/assumptions were made:
+        1. Wind speed are taken at height 100 meters above ground. However, ironmen
+           and transmission lines are typically at 50-60 meters.
+        2. Solar heat influx is set proportionally to solar short wave influx.
+        3. The incidence angle of the solar heat influx is assumed to be 90 degree.
+
+
+    Parameters
+    ----------
+    ds : xr.Dataset
+        Subset of the cutout data including all weather cells overlapping with
+        the line.
+    psi : int/float
+        Azimuth angle of the line in degree, that is the incidence angle of the line
+        with a pointer directing north (90 is east, 180 is south, 270 is west).
+    R : float
+        Resistance of the conductor in [Ω/m] at maximally allowed temperature Ts.
+    D : float,
+        Conductor diameter.
+    Ts : float
+        Maximally allowed surface temperature (typically 100°C).
+    epsilon : float
+        Conductor emissivity.
+    alpha : float
+        Conductor absorptivity.
+
+    Returns
+    -------
+    Imax
+        xr.DataArray giving the maximal current capacity per timestep in Ampere.
+
+    """
+
+    Ta = ds["temperature"]
+    Tfilm = (Ta + Ts) / 2
+    T0 = 273.15
+
+    # 1. Convective Loss, at first forced convection
+    V = ds["wnd100m"]  # typically ironmen are about 40-60 meters high
+    mu = (1.458e-6 * Tfilm ** 1.5) / (
+        Tfilm + 383.4 - T0
+    )  # Dynamic viscosity of air (13a)
+    H = ds["height"]
+    rho = (1.293 - 1.525e-4 * H + 6.379e-9 * H ** 2) / (
+        1 + 0.00367 * (Tfilm - T0)
+    )  # (14a)
+
+    reynold = D * V * rho / mu
+
+    k = (
+        2.424e-2 + 7.477e-5 * (Tfilm - T0) - 4.407e-9 * (Tfilm - T0) ** 2
+    )  # thermal conductivity
+    anglediff = ds["wnd_azimuth"] - np.deg2rad(psi)
+    Phi = np.abs(np.mod(anglediff + np.pi / 2, np.pi) - np.pi / 2)
+    K = (
+        1.194 - np.cos(Phi) + 0.194 * np.cos(2 * Phi) + 0.368 * np.sin(2 * Phi)
+    )  # wind direction factor
+
+    Tdiff = Ts - Ta
+    qcf1 = K * (1.01 + 1.347 * reynold ** 0.52) * k * Tdiff  # (3a) in [1]
+    qcf2 = K * 0.754 * reynold ** 0.6 * k * Tdiff  # (3b) in [1]
+
+    qcf = np.maximum(qcf1, qcf2)
+
+    #  natural convection
+    qcn = 3.645 * rho * D ** 0.75 * Tdiff ** 1.25
+
+    # convection loss is the max between forced and natural
+    qc = np.maximum(qcf, qcn)
+
+    # 2. Radiated Loss
+    qr = 17.8 * D * epsilon * ((Ts / 100) ** 4 - (Ta / 100) ** 4)
+
+    # 3. Solar Radiance Heat Gain
+    Q = ds["influx_direct"]  # assumption, this is short wave and not heat influx
+    A = D * 1  # projected area of conductor in square meters
+
+    qs = alpha * Q * A
+
+    Imax = np.sqrt((qc + qr - qs) / R)
+    return Imax.min("spatial") if isinstance(Imax, xr.DataArray) else Imax
+
+
+def line_rating(cutout, shapes, line_resistance, **params):
+    """
+    Create a dynamic line rating time series based on the IEEE-738 standard [1].
+
+
+    The steady-state capacity is derived from the balance between heat
+    losses due to radiation and convection, and heat gains due to solar influx
+    and conductur resistance. For more information on assumptions and modifications
+    see ``convert_line_rating``.
+
+
+    [1]“IEEE Std 738™-2012 (Revision of IEEE Std 738-2006/Incorporates IEEE Std
+        738-2012/Cor 1-2013), IEEE Standard for Calculating the Current-Temperature
+        Relationship of Bare Overhead Conductors,” p. 72.
+
+
+    Parameters
+    ----------
+    cutout : atlite.Cutout
+    shapes : geopandas.GeoSeries
+        Line shapes of the lines.
+    line_resistance : float/series
+        Resistance of the lines in Ohm/meter. Alternatively in p.u. system in
+        Ohm/1000km (see example below).
+    params : keyword arguments as float/series
+        Arguments to tweak/modify the line rating calculations based on [1].
+        Defaults are:
+            * D : 0.028 (conductor diameter)
+            * Ts : 373 (maximally allowed surface temperature)
+            * epsilon : 0.6 (conductor emissivity)
+            * alpha : 0.6 (conductor absorptivity)
+
+    Returns
+    -------
+    Current thermal limit timeseries with dimensions time x lines in Ampere.
+
+    Example
+    -------
+
+    >>> import pypsa
+    >>> import xarray as xr
+    >>> import atlite
+    >>> import numpy as np
+    >>> import geopandas as gpd
+    >>> from shapely.geometry import Point, LineString as Line
+
+    >>> n = pypsa.examples.ac_dc_meshed()
+    >>> n.calculate_dependent_values()
+    >>> x = n.buses.x
+    >>> y = n.buses.y
+    >>> buses = n.lines[["bus0", "bus1"]].values
+    >>> shapes = [Line([Point(x[b0], y[b0]), Point(x[b1], y[b1])]) for (b0, b1) in buses]
+    >>> shapes = gpd.GeoSeries(shapes, index=n.lines.index)
+
+    >>> cutout = atlite.Cutout('test', x=slice(x.min(), x.max()), y=slice(y.min(), y.max()),
+                            time='2020-01-01', module='era5', dx=1, dy=1)
+    >>> cutout.prepare()
+
+    >>> i = cutout.line_rating(shapes, n.lines.r/1e3)
+    >>> v = xr.DataArray(n.lines.v_nom, dims='name')
+    >>> s = np.sqrt(3) * i * v / 1e3 # in MW
+
+    Alternatively, the units nicely play out when we use the per unit system
+    while scaling the resistance with a factor 1000.
+
+    >>> s = np.sqrt(3) * cutout.line_rating(shapes, n.lines.r_pu * 1e3) # in MW
+
+    """
+    if not isinstance(shapes, gpd.GeoSeries):
+        shapes = gpd.GeoSeries(shapes).rename_axis("dim_0")
+
+    I = cutout.intersectionmatrix(shapes)
+    rows, cols = I.nonzero()
+
+    data = cutout.data.stack(spatial=["x", "y"])
+
+    def get_azimuth(shape):
+        coords = np.array(shape.coords)
+        start = coords[0]
+        end = coords[-1]
+        return np.arctan2(start[0] - end[0], start[1] - end[1])
+
+    azimuth = shapes.apply(get_azimuth)
+    azimuth = azimuth.where(azimuth >= 0, azimuth + np.pi)
+
+    params.setdefault("D", 0.028)
+    params.setdefault("Ts", 373)
+    params.setdefault("epsilon", 0.6)
+    params.setdefault("alpha", 0.6)
+
+    df = pd.DataFrame({"psi": azimuth, "R": line_resistance}).assign(**params)
+
+    assert df.notnull().all().all(), "Nan values encountered."
+    assert df.columns.equals(pd.Index(["psi", "R", "D", "Ts", "epsilon", "alpha"]))
+
+    dummy = xr.DataArray(np.full(len(data.time), np.nan), coords=(data.time,))
+    res = []
+    for i in range(len(df)):
+        cells_i = cols[rows == i]
+        if cells_i.size:
+            ds = data.isel(spatial=cells_i)
+            res.append(delayed(convert_line_rating)(ds, *df.iloc[i].values))
+        else:
+            res.append(dummy)
+    with ProgressBar():
+        res = compute(res)
+
+    return xr.concat(*res, dim=df.index)

atlite/cutout.py

@@ -33,7 +33,12 @@
 
 from .utils import CachedAttribute
 from .data import cutout_prepare, available_features
-from .gis import get_coords, compute_indicatormatrix, compute_availabilitymatrix
+from .gis import (
+    get_coords,
+    compute_indicatormatrix,
+    compute_availabilitymatrix,
+    compute_intersectionmatrix,
+)
 from .convert import (
     convert_and_aggregate,
     heat_demand,
@@ -44,6 +49,7 @@
     runoff,
     solar_thermal,
     soil_temperature,
+    line_rating,
 )
 from .datasets import modules as datamodules
 
@@ -527,6 +533,27 @@ def indicatormatrix(self, shapes, shapes_crs=4326):
         """
         return compute_indicatormatrix(self.grid, shapes, self.crs, shapes_crs)
 
+    def intersectionmatrix(self, shapes, shapes_crs=4326):
+        """
+        Compute the intersectionmatrix.
+
+        The intersectionmatrix is a sparse matrix with entries (i,j) being one
+        if shapes orig[j] and dest[i] are intersecting, and zero otherwise.
+
+        Note that the polygons must be in the same crs.
+
+        Parameters
+        ----------
+        orig : Collection of shapely polygons
+        dest : Collection of shapely polygons
+
+        Returns
+        -------
+        I : sp.sparse.lil_matrix
+          Intersectionmatrix
+        """
+        return compute_intersectionmatrix(self.grid, shapes, self.crs, shapes_crs)
+
     def uniform_layout(self):
         """Get a uniform capacity layout for all grid cells."""
         return xr.DataArray(1, [self.coords["y"], self.coords["x"]])
@@ -605,3 +632,5 @@ def layout_from_capacity_list(self, data, col="Capacity"):
     runoff = runoff
 
     hydro = hydro
+
+    line_rating = line_rating

atlite/datasets/era5.py

@@ -40,7 +40,7 @@ def nullcontext():
 
 features = {
     "height": ["height"],
-    "wind": ["wnd100m", "roughness"],
+    "wind": ["wnd100m", "wnd_azimuth", "roughness"],
     "influx": ["influx_toa", "influx_direct", "influx_diffuse", "albedo"],
     "temperature": ["temperature", "soil temperature"],
     "runoff": ["runoff"],
@@ -99,6 +99,9 @@ def get_data_wind(retrieval_params):
     ds["wnd100m"] = np.sqrt(ds["u100"] ** 2 + ds["v100"] ** 2).assign_attrs(
         units=ds["u100"].attrs["units"], long_name="100 metre wind speed"
     )
+    # span the whole circle: 0 is north, π/2 is east, -π is south, 3π/2 is west
+    azimuth = np.arctan2(ds["u100"], ds["v100"])
+    ds["wnd_azimuth"] = azimuth.where(azimuth >= 0, azimuth + 2 * np.pi)
     ds = ds.drop_vars(["u100", "v100"])
     ds = ds.rename({"fsr": "roughness"})
 

atlite/gis.py

@@ -153,6 +153,40 @@ def compute_indicatormatrix(orig, dest, orig_crs=4326, dest_crs=4326):
     return indicator
 
 
+def compute_intersectionmatrix(orig, dest, orig_crs=4326, dest_crs=4326):
+    """
+    Compute the intersectionmatrix.
+
+    The intersectionmatrix is a sparse matrix with entries (i,j) being one
+    if shapes orig[j] and dest[i] are intersecting, and zero otherwise.
+
+    Note that the polygons must be in the same crs.
+
+    Parameters
+    ----------
+    orig : Collection of shapely polygons
+    dest : Collection of shapely polygons
+
+    Returns
+    -------
+    I : sp.sparse.lil_matrix
+      Intersectionmatrix
+    """
+    orig = orig.geometry if isinstance(orig, gpd.GeoDataFrame) else orig
+    dest = dest.geometry if isinstance(dest, gpd.GeoDataFrame) else dest
+    dest = reproject_shapes(dest, dest_crs, orig_crs)
+    intersection = sp.sparse.lil_matrix((len(dest), len(orig)), dtype=float)
+    tree = STRtree(orig)
+    idx = dict((id(o), i) for i, o in enumerate(orig))
+
+    for i, d in enumerate(dest):
+        for o in tree.query(d):
+            j = idx[id(o)]
+            intersection[i, j] = o.intersects(d)
+
+    return intersection
+
+
 class ExclusionContainer:
     """Container for exclusion objects and meta data."""
 

atlite/wind.py

@@ -9,6 +9,7 @@
 """
 
 import numpy as np
+import re
 
 import logging
 
@@ -59,7 +60,7 @@ def extrapolate_wind_speed(ds, to_height, from_height=None):
 
     if from_height is None:
         # Determine closest height to to_name
-        heights = np.asarray([int(s[3:-1]) for s in ds if s.startswith("wnd")])
+        heights = np.asarray([int(s[3:-1]) for s in ds if re.match(r"wnd\d+m", s)])
 
         if len(heights) == 0:
             raise AssertionError("Wind speed is not in dataset")