examples/working-with-csp.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "431006d2-8530-496d-b82c-3dd979fcfdb3",
   "metadata": {},
   "source": [
    "# Concentrated Solar Power in atlite\n",
    "\n",
    "## Introduction\n",
    "\n",
    "In this page we will have a look at how `atlite` can be used to determine the output from CSP (Concentrated Solar Power).\n",
    "`atlite` currently implements two CSP technologies:\n",
    "\n",
    "* Parabolic trough\n",
    "* Solar tower\n",
    "\n",
    "For both technologies `atlite` can be used to determine the heat output from the heliostat field or solar receiver.\n",
    "Subsequent steps, e.g. heat storage, electricity conversion and dispatch are not covered by atlite and have to be handled\n",
    "by an external model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "37d30d6b-0b2a-4736-9bfc-8315f54f69fc",
   "metadata": {},
   "outputs": [],
   "source": [
    "import geopandas as gpd\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import xarray as xr\n",
    "import yaml\n",
    "\n",
    "import atlite"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5052eb8f-5ffb-4249-97f3-7dca75b7fd77",
   "metadata": {},
   "source": [
    "### Technologies\n",
    "The two technologies are different in the way that they absorb and convert different parts of the direct radiation hitting the surface.\n",
    "\n",
    "* **Parabolic troughs** only have horizontal tracking of the solar position and thus is modelled to only absorb the **Direct Horizontal Irradiation**, i.e. direct radiation hitting the horizontal plane.\n",
    "* **Solar tower** heliostats on the other hand require two-axis tracking of the solar position in order to reflect the irradiation into the absorber; this technology is therefore modelled to  absorb the **Direct Normal Irradiation**."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ffe2a0c0-1baf-4194-bf52-21edb17e772c",
   "metadata": {},
   "source": [
    "### Installation specific efficiencies and properties\n",
    "\n",
    "`atlite` uses a simplified approach in modelling CSP. The technology-specific direct irradiation is modified by an efficiency which depends on the local solar position.\n",
    "The reason for an efficiency which is solar dependend is that a CSP installation can not perfectly absorb all direct irradiation at any given time.\n",
    "\n",
    "Effects causing less-than-perfect efficiency are e.g.:\n",
    "\n",
    "* shadowing between rows and adjacents heliostats\n",
    "* maximum tracking angles\n",
    "* system design for reduced efficiency during peak irradiation hours in favour of increased efficiency during off-peak sunset/sunrise\n",
    "\n",
    "The efficiency can be very specific to an CSP installation and its design.\n",
    "Therefore `atlite` uses the properties of specific CSP installations - similar as for wind turbines or PV panels.\n",
    "`atlite` ships with these three reference CSP installations:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "189304ad-8890-4a84-a311-56e64ece2ec0",
   "metadata": {},
   "outputs": [],
   "source": [
    "atlite.cspinstallations.keys()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "27765b99-7e69-4479-ad1a-08dc21056134",
   "metadata": {},
   "source": [
    "An installation is characterised by a `dict` of properties"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "433e917a-26fe-4705-a93f-211e0c0e8bf7",
   "metadata": {},
   "outputs": [],
   "source": [
    "atlite.resource.get_cspinstallationconfig(atlite.cspinstallations.SAM_parabolic_trough)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0c5e1859-1002-49b1-9dba-ed95d92baf14",
   "metadata": {},
   "source": [
    "\n",
    "Explanation of keys and values:\n",
    "\n",
    "| key | value | description |\n",
    "| --- | ----- | ----------- |\n",
    "| name | str | The installations name. |\n",
    "| source | str | Source of the installation details. |\n",
    "| [technology](#technologies) | One of {'parabolic trough', 'solar tower', None}. | Technology to use for this installation. If 'None', technology has to be specified later. |\n",
    "| r_irradiance | number | Reference irradiance for this installation. Used to calculated capacity factor and maximum output of system. Unit: W/m^2 .|\n",
    "| [efficiency](#efficiency) | xarray.DataArray | DataArray containing solar position dependend conversion efficiency between Normal Irradiation and thermal output. |\n",
    "| path | Path | (automatically added) Path the information was loaded from. |\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bf9d8aa2-4aeb-40ca-bc17-4d3675ddf251",
   "metadata": {
    "tags": []
   },
   "source": [
    "### Efficiency"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f1c2a03c-ba23-4068-81b6-ce54bdc3f475",
   "metadata": {},
   "source": [
    "The efficiencies behaviour of the three installations are quite different"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7d17418f-c177-4033-b219-7ffa9f2cca73",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(1, 3, figsize=(20, 5))\n",
    "\n",
    "for i, name in enumerate(atlite.cspinstallations.keys()):\n",
    "    config = atlite.resource.get_cspinstallationconfig(atlite.cspinstallations[name])\n",
    "\n",
    "    config[\"efficiency\"].plot(\n",
    "        x=\"azimuth [deg]\",\n",
    "        y=\"altitude [deg]\",\n",
    "        cmap=\"plasma\",\n",
    "        ax=axes[i],\n",
    "        vmin=0,\n",
    "        vmax=1.0,\n",
    "    )\n",
    "    axes[i].set_title(f\"{name} (technology: {config['technology']})\")\n",
    "    axes[i].set_xlim(0, 360)\n",
    "    axes[i].set_ylim(0, 90)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4926b097-7875-4f00-addc-93fa39256b8e",
   "metadata": {},
   "source": [
    "For the two 'SAM_*' models the efficiency strongly depends on the solar azimuth and altitude.\n",
    "Due to the installations' specifics, no efficiency is defined for high solar azimuth positions.\n",
    "These installations may therefore not be suited for locations with high solar azimuths.\n",
    "Both models also have technologies predefined, such that the 'technology' parameter does not need to be supplied during handling.\n",
    "The installation details for 'SAM_parabolic_trough' and 'SAM_solar_tower' were retrieved and modelled with NREL's System Advisor Model.\n",
    "For details on the process see [the section below](#extracting-efficiencies-from-sam).\n",
    "\n",
    "The 'lossless_installation' is an universal installation which has by default perfect efficiency and is thus able to convert the direct irradiation of either technology into an heat output.\n",
    "It can be used to easily create fixed-efficiency installations, e.g. an installation with solar to heat efficiency of 65% may be created as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f3568c86-7e23-456c-8944-886ad0beb41f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load installation configuration by\n",
    "config = atlite.resource.get_cspinstallationconfig(\"lossless_installation\")\n",
    "# or\n",
    "config = atlite.resource.get_cspinstallationconfig(\n",
    "    atlite.cspinstallations.lossless_installation\n",
    ")\n",
    "\n",
    "# Reduce efficiency to 65%\n",
    "config[\"efficiency\"] *= 0.65"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "215e9e97-b93f-473a-bab4-b92b62f2abc0",
   "metadata": {},
   "source": [
    "## Resource evaluation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "98458c18-46ff-4477-b474-a6968b8776ca",
   "metadata": {},
   "source": [
    "To evaluate the resources and suitability of locations for CSP, we first need to [create an cutout](https://atlite.readthedocs.io/en/latest/examples/create_cutout.html) for the area of interest:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cd370f1c-7cc0-4cb7-952b-ca5d56dd4924",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "world = gpd.read_file(gpd.datasets.get_path(\"naturalearth_lowres\"))\n",
    "spain = world.query('name == \"Spain\"')\n",
    "\n",
    "cutout = atlite.Cutout(\n",
    "    path=\"spain_2017_era5.nc\",\n",
    "    module=\"era5\",\n",
    "    bounds=spain.iloc[0].geometry.bounds,\n",
    "    time=\"2017\",\n",
    ")\n",
    "\n",
    "cutout.prepare([\"influx\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b044579a-ed42-42ed-840b-12a5ed5099b6",
   "metadata": {},
   "source": [
    "### Capacity factor"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dbddb6b8",
   "metadata": {},
   "source": [
    "We can then calculate the capacity factor for each grid cell by specifying an installation and (if not included in the installation configuration) a 'technology':"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bafabada-e683-4fc9-b281-72f999e88960",
   "metadata": {},
   "outputs": [],
   "source": [
    "cf = cutout.csp(installation=\"SAM_solar_tower\", capacity_factor=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "05a13da5-70f7-4bbd-9497-6f4e2f43dbc0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot capacity factor with Spain shape overlay\n",
    "fig, ax = plt.subplots()\n",
    "cf.plot(cmap=\"inferno\", ax=ax)\n",
    "spain.plot(fc=\"None\", ec=\"black\", ax=ax)\n",
    "fig.suptitle(\"CSP capacity factor\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "123230d8-beaf-42e8-9fda-74f03cbfbe9c",
   "metadata": {},
   "source": [
    "Alternatively the specific generation can be calculated by setting `capacity_factor=False` (default).\n",
    "An comparison between the `SAM_solar_tower` installation and `lossless_installation` shows the difference due to the solar position dependend efficiency:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "608508da-5b90-41b6-955c-ac417f186ab9",
   "metadata": {},
   "source": [
    "We take as comparison an existing installation of a Solar Tower, the [PS10 in Spain](https://en.wikipedia.org/wiki/PS10_solar_power_plant)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c374d216-d30d-4d6f-8ea6-3fcad533fea7",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Calculation of Capacity Factor and Specific Generation for: SAM_solar_tower installation\n",
    "st = {\n",
    "    \"capacity factor\": cutout.csp(\n",
    "        installation=\"SAM_solar_tower\", capacity_factor=True\n",
    "    ).rename(\"SAM_solar_tower CF\"),\n",
    "    \"specific generation\": cutout.csp(installation=\"SAM_solar_tower\").rename(\n",
    "        \"SAM_solar_tower SG\"\n",
    "    ),\n",
    "}\n",
    "\n",
    "# Calculation of Capacity Factor and Specific Generation for: lossless solar tower installation\n",
    "ll = {\n",
    "    \"capacity factor\": cutout.csp(\n",
    "        installation=\"lossless_installation\",\n",
    "        technology=\"solar tower\",\n",
    "        capacity_factor=True,\n",
    "    ).rename(\"lossless_installation CF\"),\n",
    "    \"specific generation\": cutout.csp(\n",
    "        installation=\"lossless_installation\",\n",
    "        technology=\"solar tower\",\n",
    "    ).rename(\"lossless_installation SG\"),\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0c9de37e-8f80-4020-9502-dd33d05b005d",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Plot results side by side\n",
    "fig, axes = plt.subplots(2, 2, figsize=(16, 10))\n",
    "\n",
    "st[\"capacity factor\"].plot(ax=axes[0][0], cmap=\"inferno\")\n",
    "st[\"specific generation\"].plot(ax=axes[1][0], cmap=\"inferno\")\n",
    "axes[0][0].set_title(\"SAM_solar_tower\")\n",
    "\n",
    "ll[\"capacity factor\"].plot(ax=axes[0][1], cmap=\"inferno\")\n",
    "ll[\"specific generation\"].plot(ax=axes[1][1], cmap=\"inferno\")\n",
    "axes[0][1].set_title(\"lossless_installation\")\n",
    "\n",
    "# Overlay Spainish borders\n",
    "for ax in axes.ravel():\n",
    "    spain.plot(ax=ax, fc=\"none\", ec=\"black\")\n",
    "\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "368bcef2-73a2-4bf1-aa71-aa313e4505e1",
   "metadata": {},
   "source": [
    "### Capacity layout for a specific location"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "af092fac-2eb1-4028-8ebf-2b1bd79e5aa5",
   "metadata": {},
   "source": [
    "We take as comparison an existing installation of a Solar Tower, the [PS10 in Spain](https://en.wikipedia.org/wiki/PS10_solar_power_plant):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1a64a6bd-86a6-4d9d-8ac8-da703c1b325b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Layout with different installed capacity (due to different SAM model)\n",
    "config = atlite.resource.get_cspinstallationconfig(\"SAM_solar_tower\")\n",
    "# solar field: 624 mirrors with 120m^2 each\n",
    "area = 624 * 120  # 74880 m^2\n",
    "\n",
    "# installed power = 950 W/m^2 * area = 90.2 MW\n",
    "installed_power = config[\"r_irradiance\"] * area / 1.0e6\n",
    "\n",
    "\n",
    "# actual PS10 location\n",
    "# see https://geohack.toolforge.org/geohack.php?pagename=PS10_solar_power_plant&params=37_26_32_N_06_15_15_W_type:landmark_region:ES\n",
    "location = {\"x\": -6.254167, \"y\": 37.442222}\n",
    "\n",
    "# Determine location on cutout grid\n",
    "nearest_location = {\n",
    "    v: cutout.grid[v].iloc[cutout.grid[v].sub(location[v]).abs().idxmin()]\n",
    "    for v in [\"x\", \"y\"]\n",
    "}\n",
    "\n",
    "layout = xr.zeros_like(cf)\n",
    "layout.loc[dict(x=nearest_location[\"x\"], y=nearest_location[\"y\"])] = installed_power"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "17cdadd9-3e9f-4d04-b0ec-c958e01f1627",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot layout with  Spain shape overlay\n",
    "fig, ax = plt.subplots()\n",
    "layout.plot(cmap=\"inferno_r\", ax=ax)\n",
    "spain.plot(fc=\"None\", ec=\"black\", ax=ax)\n",
    "fig.suptitle(\"Layout with PS10\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ade70956-db16-4794-bee1-ee4eacc0d321",
   "metadata": {},
   "source": [
    "### Solar tower: Hourly time-series and monthly output"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8ea13308-5c12-42ad-b94d-265a18d3bfff",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate time-series for layout with both installation configurations\n",
    "time_series = xr.merge(\n",
    "    [\n",
    "        cutout.csp(\n",
    "            installation=\"lossless_installation\",\n",
    "            technology=\"solar tower\",\n",
    "            layout=layout,\n",
    "        ).rename(\"lossless_installation\"),\n",
    "        cutout.csp(installation=\"SAM_solar_tower\", layout=layout).rename(\n",
    "            \"SAM_solar_tower\"\n",
    "        ),\n",
    "    ]\n",
    ")\n",
    "\n",
    "# Load reference time-series from file\n",
    "df = pd.read_csv(\"../profiles_and_efficiencies_from_sam/ST-salt_time-series_spain.csv\")\n",
    "df[\"time\"] = time_series[\"time\"].values + pd.Timedelta(\n",
    "    1, unit=\"h\"\n",
    ")  # account for reference values in UTC-1 instead of UTC\n",
    "df = df.set_index(\"time\")\n",
    "df = df[\"Estimated receiver thermal power TO HTF | (MWt)\"] * (\n",
    "    # Rescale: Simulation in SAM uses solar field of 1 269 054 m^2 and 950 W/m^2 reference irradiation = 1 205.6 MW\n",
    "    installed_power / (950 * 1269054 / 1e6)\n",
    ")\n",
    "time_series[\"NREL SAM\"] = df.to_xarray()\n",
    "\n",
    "time_series = time_series.squeeze().drop(\"dim_0\")  # Remove excess dimension"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "afbb6444-0d57-4a77-ac08-6c27424701c8",
   "metadata": {},
   "outputs": [],
   "source": [
    "time_series.to_array().sel(time=slice(\"2017-03-01\", \"2017-03-07\")).plot.line(\n",
    "    x=\"time\", hue=\"variable\"\n",
    ")\n",
    "\n",
    "plt.ylabel(\"Thermal output [MWh_th]\")\n",
    "plt.title(\"Example output time-series\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cd2c4b57-f790-4502-8119-fb3e5235742a",
   "metadata": {},
   "source": [
    "From the figure above the difference between perfect efficiency ('lossless_installation') and solar position dependend conversion efficiency ('SAM_solar_tower') can be seen.\n",
    "Mostly for low solar altitudes (where the DNI calculation may overestimate DNI due to cosine errors), the imperfect efficiency shows a smoother generation profile than the perfect installation.\n",
    "\n",
    "The 'reference' data included in the figure shows the output of the more sophisticated SAM CSP model for the same location and installation.\n",
    "In general a qualitatively acceptable fit is achieved with the simplified approach used in `atlite`.\n",
    "\n",
    "Looking at the monthly energy output we see an approximate match between the `atlite` and reference results:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cdbadb45-e867-4fc9-b6b7-2021391e5a59",
   "metadata": {},
   "outputs": [],
   "source": [
    "time_series.groupby(\"time.month\").sum().to_dataframe().plot.bar()\n",
    "plt.ylabel(\"Thermal output [MWh_th]\")\n",
    "plt.title(\"Monthly output\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7ab3230b-ab3d-4672-9b2a-14d0cf6aad91",
   "metadata": {},
   "source": [
    "### Differences between NREL's SAM and atlite output"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dca900e2-c0f1-4b64-9dbc-1dda8f5a305e",
   "metadata": {},
   "source": [
    "Differences between the `atlite` results and more sophisticated SAM reference model can be due to a multitude of effects, e.g.:\n",
    "\n",
    "* The reference model uses weather data from the National Solar Radiation Database (NSRDB), the simulated results from atlite use ERA5 data. NSRDB relies on a physical model which combines satellite data with reanalysis data, compared to ERA5 which is a pure reanalysis dataset.\n",
    "* `atlite` uses an simple approach to determining DNI for solar tower heliostat fields\n",
    "* The `reference` data against which the `atlite` output is compared against reflects the thermal HTF (heat transfer fluid) output; the `atlite` numbers on the other hand reflect the absorbed thermal energy\n",
    "* The reference model accounts for wind velocities and necessary heliostat stowing\n",
    "* The reference model takes maintenance and cleaning outages into account\n",
    "* The reference model can reduce heliostat performance due to dispatch reasons"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "770e1482-2f4f-411e-a4c4-75969ccde376",
   "metadata": {},
   "source": [
    "### Parabolic trough: Hourly time-series and monthly output"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bccb09fa-12f4-441d-bff1-f5f8a6644ebe",
   "metadata": {},
   "source": [
    "The issues of overestimating direct irradiation does not occur with `technology=parabolic trough`, in which case an DNI esimation is not necessary and `atlite` can directly rely on DHI retrieved from the cutout's dataset.\n",
    "The following calculation and figure show this using a fictive parabolic trough CSP plant in the same location as the PS10 plant from above:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cf28821f-0ab3-4918-9b2c-2c7ee67c6ac3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Layout with different installed capacity (due to different SAM model)\n",
    "config = atlite.resource.get_cspinstallationconfig(\"SAM_parabolic_trough\")\n",
    "# solar field size in m^2 of fictive plant\n",
    "area = 881664\n",
    "\n",
    "# installed power = 950 W/m^2 * area = 1205.0 MW\n",
    "installed_power = config[\"r_irradiance\"] * area / 1.0e6\n",
    "layout.loc[dict(x=nearest_location[\"x\"], y=nearest_location[\"y\"])] = installed_power\n",
    "\n",
    "# Calculate time-series for layout with both installation configurations\n",
    "time_series = xr.merge(\n",
    "    [\n",
    "        cutout.csp(\n",
    "            installation=\"lossless_installation\",\n",
    "            technology=\"parabolic trough\",\n",
    "            layout=layout,\n",
    "        ).rename(\"lossless_installation\"),\n",
    "        cutout.csp(installation=\"SAM_parabolic_trough\", layout=layout).rename(\n",
    "            \"SAM_parabolic_trough\"\n",
    "        ),\n",
    "    ]\n",
    ")\n",
    "\n",
    "# Load reference time-series from file\n",
    "df = pd.read_csv(\n",
    "    \"../profiles_and_efficiencies_from_sam/PT-physical_time-series_spain.csv\"\n",
    ")\n",
    "# account for reference values in UTC-1 instead of UTC\n",
    "df[\"time\"] = time_series[\"time\"].values + pd.Timedelta(1, unit=\"h\")\n",
    "\n",
    "df = df.set_index(\"time\")\n",
    "df = df[\"Estimate rec. thermal power TO HTF | (MWt)\"]\n",
    "time_series[\"NREL SAM\"] = df.to_xarray()\n",
    "\n",
    "time_series = time_series.squeeze().drop(\"dim_0\")  # Remove excess dimension"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c7c47f73-1cc4-4793-a34c-df62145e8495",
   "metadata": {},
   "outputs": [],
   "source": [
    "time_series.to_array().sel(time=slice(\"2017-03-01\", \"2017-03-07\")).plot.line(\n",
    "    x=\"time\", hue=\"variable\"\n",
    ")\n",
    "plt.ylabel(\"Thermal output [MWh_th]\")\n",
    "plt.title(\"Example output time-series\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aec54683",
   "metadata": {},
   "source": [
    "For `technology=\"parabolic trough\"`  the monthly output is also better:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "029df2da-c581-400b-9d42-903fb2f2095d",
   "metadata": {},
   "outputs": [],
   "source": [
    "time_series.groupby(\"time.month\").sum().to_dataframe().plot.bar()\n",
    "plt.ylabel(\"Thermal output [MWh_th]\")\n",
    "plt.title(\"Monthly output\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4d721c1d-0cdc-46a2-a3c5-734a8866c7f9",
   "metadata": {},
   "source": [
    "## ERA5 vs. SARAH2 data"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5c189c17-6763-4201-947b-d05f39a03a35",
   "metadata": {},
   "source": [
    "Next to ERA5 data, SARAH2 data can also be used as the underpinning `cutout` data.\n",
    "SARAH2 is generally considered more accurate for radiation products as it uses satellite data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7eac0412-eee8-4017-98ab-bc6ed740b96f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create cutout using SARAH data\n",
    "cutout_sarah = atlite.Cutout(\n",
    "    path=\"spain_2017_sarah.nc\",\n",
    "    module=\"sarah\",\n",
    "    bounds=spain.iloc[0].geometry.bounds,\n",
    "    time=\"2017\",\n",
    ")\n",
    "\n",
    "# Save space and computing time: only need the influx feature for demonstration\n",
    "cutout_sarah.prepare(features=[\"influx\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bcc0f41f-c0c5-4132-860c-8b986550063a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate time-series for layout with both installation configurations\n",
    "time_series = xr.merge(\n",
    "    [\n",
    "        cutout_sarah.csp(installation=\"SAM_parabolic_trough\", layout=layout).rename(\n",
    "            \"SARAH\"\n",
    "        ),\n",
    "        cutout.csp(\n",
    "            installation=\"SAM_parabolic_trough\",\n",
    "            layout=layout,\n",
    "        ).rename(\"ERA5\"),\n",
    "    ]\n",
    ")\n",
    "\n",
    "# Load reference NREL SAM time-series from file\n",
    "df = pd.read_csv(\n",
    "    \"../profiles_and_efficiencies_from_sam/PT-physical_time-series_spain.csv\"\n",
    ")\n",
    "\n",
    "# account for reference values in UTC-1 instead of UTC\n",
    "df[\"time\"] = time_series[\"time\"].values + pd.Timedelta(1, unit=\"h\")\n",
    "\n",
    "df = df.set_index(\"time\")\n",
    "df = df[\"Estimate rec. thermal power TO HTF | (MWt)\"]\n",
    "time_series[\"NREL SAM\"] = df.to_xarray()\n",
    "\n",
    "time_series = time_series.squeeze().drop(\"dim_0\")  # Remove excess dimension"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "648a8d17-d376-4e10-adba-68af63df3da4",
   "metadata": {},
   "outputs": [],
   "source": [
    "time_series.to_array().sel(time=slice(\"2017-03-01\", \"2017-03-07\")).plot.line(\n",
    "    x=\"time\", hue=\"variable\"\n",
    ")\n",
    "plt.ylabel(\"Thermal output [MWh_th]\")\n",
    "plt.title(\"Example output time-series\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "44f5e7c1-3495-4e13-a8c4-d6df70c968f8",
   "metadata": {},
   "source": [
    "In the monthly aggregated output it can be seen that the SARAH dataset has a lower overestimation error compared to the ERA5 dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1d546103-3ff6-4420-9f00-098ae6473ec2",
   "metadata": {},
   "outputs": [],
   "source": [
    "time_series.groupby(\"time.month\").sum().to_dataframe().plot.bar()\n",
    "plt.ylabel(\"Thermal output [MWh_th]\")\n",
    "plt.title(\"Monthly output\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "14580e09-2975-4c99-935f-689cc9f7678f",
   "metadata": {},
   "source": [
    "Neither dataset (ERA5 and SARAH) are replicating the generation output of a singular, specific plant modelled with NREL's SAM using NSRDB data.\n",
    "`atlite` generally is not designed with modelling generation output of singular plants.\n",
    "Rather `atlite` is meant to be used to assess the ressources of a larger geographic extent for potential CSP output.\n",
    "\n",
    "We end this comparison with a sometimes used, very simple approach of only using the DNI time-series to model the CSP generation.\n",
    "Compared on a monthly basis the `atlite` output is closer to the NREL SAM reference model.\n",
    "For time-series output one can compare the time-series of the \"lossless\" plant configuration shown earlier.\n",
    "Generally, the qualitative modelling performance of this approach is poorer compared to the `atlite` or `SAM` methods."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "63899423-154c-462c-ae56-14b5cf19a83f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Monthly average DNI data for a TMY from GlobalSolarAtlas\n",
    "# URL: https://globalsolaratlas.info/detail?s=37.5,-6.3&m=site&c=33.583674,-10.898438,3 (accessed: 2021-11-22)\n",
    "df = pd.DataFrame(\n",
    "    [122.5, 131, 161.9, 169.1, 197.1, 230.3, 259.6, 232.8, 175.3, 142.7, 122.8, 116.6],\n",
    "    index=pd.Index(np.arange(1, 13, 1), name=\"month\"),\n",
    "    columns=[\"GlobalSolarAtlas\"],\n",
    ")\n",
    "\n",
    "ds = time_series.groupby(\"time.month\").sum() / area * 1e3  # in kWh/m^2\n",
    "\n",
    "ds = ds.to_dataframe()\n",
    "\n",
    "ds[\"GlobalSolarAtlas\"] = df[\"GlobalSolarAtlas\"]\n",
    "\n",
    "ds.plot.bar()\n",
    "plt.ylabel(\"Thermal output [kWh/m^2 heliostat area]\")\n",
    "plt.title(\"Monthly output\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8e26003a-632d-4a2c-b716-d77625dcce67",
   "metadata": {},
   "source": [
    "## Extracting efficiencies from SAM"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c4ee173f-9832-44b3-8afb-2d2f591a85f5",
   "metadata": {},
   "source": [
    "Instead of implementing a detailed physical model for CSP plants in `atlite` an approach was selected where an efficiency map is handed over to `atlite` for simulating the conversion efficiency and processes inside the solar field of an CSP plant.\n",
    "The efficiency maps which come with `atlite` were created using [NREL's System Advisor Model](https://sam.nrel.gov/), which allows for a detailed specification and simulation of CSP plants (solar tower and parabolic trough).\n",
    "\n",
    "The hourly solar field efficiencies and solar positions were then exported to a `.csv` file.\n",
    "Using this import, the efficiency maps for `atlite` were then created.\n",
    "In the following these steps are shown and the output is compared to the input from SAM."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "65a41892-f94a-4bc7-b051-4b45ca042019",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Read SAM output: Hourly solar field efficiency and solar positions\n",
    "df = pd.read_csv(\"../profiles_and_efficiencies_from_sam/PT-physical.csv\")\n",
    "\n",
    "df = df.rename(\n",
    "    columns={\n",
    "        \"Hourly Data: Resource Solar Azimuth (deg)\": \"azimuth\",\n",
    "        \"Hourly Data: Resource Solar Zenith (deg)\": \"zenith (deg)\",\n",
    "        \"Hourly Data: Field optical efficiency before defocus\": \"efficiency\",\n",
    "    }\n",
    ")\n",
    "\n",
    "# Convert solar zenith to solar altitude\n",
    "df[\"altitude\"] = 90 - df[\"zenith (deg)\"]\n",
    "\n",
    "# Only non-zero entries from hourly data are relevant\n",
    "df = df[df[\"efficiency\"] > 0.0]\n",
    "\n",
    "# Efficiency in %, avoids dropping significant places in round(1) later\n",
    "df[\"efficiency\"] *= 100\n",
    "\n",
    "# Reduce noise in solar position by rounding each position to 2°\n",
    "# Creates a lot of duplicates\n",
    "base = 2\n",
    "\n",
    "df[[\"azimuth\", \"altitude\"]] = (df[[\"azimuth\", \"altitude\"]] / base).round(0) * base\n",
    "\n",
    "# Drop duplicates: Keep only highest efficiency per position\n",
    "df = df.groupby([\"altitude\", \"azimuth\"]).max()\n",
    "\n",
    "# Drop all columns except for efficiency\n",
    "df = df[[\"efficiency\"]]\n",
    "\n",
    "da = df.to_xarray()[\"efficiency\"]\n",
    "\n",
    "# Interpolate values to a finer grid and fill missing values by extrapolation\n",
    "# Order is relevant: Start with Azimuth (where we have sufficient values) and then continue with altitude\n",
    "da = (\n",
    "    da.interpolate_na(\"azimuth\")\n",
    "    .interpolate_na(\"altitude\")\n",
    "    .interpolate_na(\"azimuth\", fill_value=\"extrapolate\")\n",
    ")\n",
    "\n",
    "# Use rolling horizon to smooth values, average over 3x3 adjacent values per pixel\n",
    "da = (\n",
    "    da.rolling(azimuth=3, altitude=3)\n",
    "    .mean()\n",
    "    .interpolate_na(\"altitude\", fill_value=\"extrapolate\")\n",
    "    .interpolate_na(\"azimuth\", fill_value=\"extrapolate\")\n",
    ")\n",
    "\n",
    "# Create second dataset, mirrored around 90° (covering 90° to -90° with same numbers)\n",
    "da = da.sel(azimuth=slice(90, 270))\n",
    "da_new = da.assign_coords(azimuth=(((da.azimuth - 180) * (-1)) % 360)).sortby(\"azimuth\")\n",
    "da = xr.merge([da, da_new])[\"efficiency\"]\n",
    "\n",
    "# Reduce resolution to coarser grid for reduced storage footprint\n",
    "interpolation_step = 5\n",
    "da = da.interp(\n",
    "    azimuth=np.arange(da[\"azimuth\"].min(), da[\"azimuth\"].max(), interpolation_step),\n",
    "    altitude=np.arange(0, 91, interpolation_step),\n",
    "    kwargs={\"fill_value\": \"extrapolate\"},\n",
    ").clip(min=0.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b565c947-598a-41d3-9fd5-4cde4195579a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot efficiencies against each other\n",
    "fig, axes = plt.subplots(ncols=3, nrows=1, sharey=True, figsize=(30, 5))\n",
    "\n",
    "ax = axes[0]\n",
    "df.to_xarray()[\"efficiency\"].plot(cmap=\"viridis\", ax=ax)\n",
    "ax.set_title(\"Imported from SAM\")\n",
    "ax.set_xlim(0, 360)\n",
    "\n",
    "ax = axes[1]\n",
    "da.plot(cmap=\"viridis\", ax=ax)\n",
    "ax.set_title(\"Efficiencies interpolated+smoothed for atlite\")\n",
    "ax.set_xlim(0, 360)\n",
    "\n",
    "ax = axes[2]\n",
    "(da.interp_like(df.to_xarray()) - df.to_xarray()[\"efficiency\"]).plot(\n",
    "    cmap=\"viridis\", ax=ax\n",
    ")\n",
    "ax.set_title(\"Absolute difference between SAM and atlite\")\n",
    "ax.set_xlim(0, 360)\n",
    "\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3bf1db17-571b-4cd5-9211-57cf19c13f40",
   "metadata": {},
   "source": [
    "The differences between the original efficiency map from SAM and the derived efficiency map for atlite are low, in most cases below 1%.\n",
    "The efficiency map is then written to a `.yaml` file, combined with the remaining installation information for the specific CSP plant (also taken from SAM)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aeb9076e-1b6c-4dac-bcd6-a478c6350b81",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = da.to_dataframe().reset_index()\n",
    "df = df.rename(columns={\"efficiency\": \"value\"})\n",
    "df[[\"altitude\", \"azimuth\"]] = df[[\"altitude\", \"azimuth\"]].astype(int)\n",
    "df = df.to_dict(\"dict\")\n",
    "\n",
    "config = {\n",
    "    \"name\": \"Parabolic Trough field\",\n",
    "    \"source\": \"NREL System Advisor Model\",\n",
    "    \"technology\": \"parabolic trough\",\n",
    "    \"r_irradiance\": 950,  # W/m2,\n",
    "    \"efficiency\": df,\n",
    "}\n",
    "\n",
    "with open(\"SAM_solar_tower.yaml\", \"w\") as f:\n",
    "    yaml.safe_dump(config, f, default_flow_style=False, sort_keys=False)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "",
   "language": "python",
   "name": ""
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  },
  "toc-autonumbering": false,
  "toc-showcode": false,
  "toc-showmarkdowntxt": false,
  "toc-showtags": false
 },
 "nbformat": 4,
 "nbformat_minor": 5
}