Add land cover analysis for MD, UA (Quick & Dirty).

Snakefile

@@ -66,6 +66,43 @@ if config['enable'].get('retrieve_databundle', True):
         log: "logs/retrieve_databundle.log"
         script: 'scripts/retrieve_databundle.py'
 
+# Downloading Copernicus Global Land Cover for land cover and land use:
+# Website: https://land.copernicus.eu/global/products/lc
+rule download_copernicus_land_cover:
+    input:
+        HTTP.remote(
+            "zenodo.org/record/3939050/files/PROBAV_LC100_global_v3.0.1_2019-nrt_Discrete-Classification-map_EPSG-4326.tif",
+            static=True,
+        ),
+    output:
+        "resources/Copernicus_LC100_global_v3.0.1_2019-nrt_Discrete-Classification-map_EPSG-4326.tif",
+    shell:
+        "mv {input} {output}"
+
+rule determine_availability_matrix_MD_UA:
+    input:
+        copernicus="resources/Copernicus_LC100_global_v3.0.1_2019-nrt_Discrete-Classification-map_EPSG-4326.tif",
+        gebco=lambda w: ("data/bundle/GEBCO_2014_2D.nc"
+                         if "max_depth" in config["renewable"][w.technology].keys()
+                         else []),
+        country_shapes='resources/country_shapes.geojson',
+        offshore_shapes='resources/offshore_shapes.geojson',
+        regions=lambda w: ("resources/regions_onshore.geojson"
+                                   if w.technology in ('onwind', 'solar')
+                                   else "resources/regions_offshore.geojson"),
+        cutout=lambda w: "cutouts/" + config["renewable"][w.technology]['cutout'] + ".nc"
+    output:
+        availability_matrix="resources/availability_matrix_MD-UA_{technology}.nc",
+    log:
+        "logs/determine_availability_matrix_MD_UA_{technology}.log",
+    benchmark:
+        "benchmarks/determine_availability_matrix_MD_UA_{technology}.log",
+    threads:
+        ATLITE_NPROCESSES
+    resources:
+        mem_mb=ATLITE_NPROCESSES * 5000
+    notebook:
+        "scripts/determine_availability_matrix_MD_UA.py.ipynb"
 
 rule retrieve_load_data:
     input: HTTP.remote("data.open-power-system-data.org/time_series/2019-06-05/time_series_60min_singleindex.csv", keep_local=True, static=True)
@@ -182,6 +219,14 @@ if config['enable'].get('retrieve_natura_raster', True):
         run: move(input[0], output[0])
 
 
+# Optional input when having Ukraine (UA) or Moldova (MD) in the countries list
+if {"UA", "MD"}.intersection(set(config["countries"])):
+    opt = {
+        "availability_matrix_MD_UA":"resources/availability_matrix_MD-UA_{technology}.nc"
+    }
+else:
+    opt = {}
+
 rule build_renewable_profiles:
     input:
         base_network="networks/base.nc",
@@ -195,7 +240,8 @@ rule build_renewable_profiles:
         regions=lambda w: ("resources/regions_onshore.geojson"
                                    if w.technology in ('onwind', 'solar')
                                    else "resources/regions_offshore.geojson"),
-        cutout=lambda w: "cutouts/" + config["renewable"][w.technology]['cutout'] + ".nc"
+        cutout=lambda w: "cutouts/" + config["renewable"][w.technology]['cutout'] + ".nc",
+        **opt,
     output: profile="resources/profile_{technology}.nc",
     log: "logs/build_renewable_profile_{technology}.log"
     benchmark: "benchmarks/build_renewable_profiles_{technology}"

scripts/build_renewable_profiles.py

@@ -260,6 +260,12 @@ if __name__ == '__main__':
     else:
         availability = cutout.availabilitymatrix(regions, excluder, **kwargs)
 
+    # For Moldova and Ukraine: Overwrite parts not covered by Corine with
+    # externally determined available areas
+    if "availability_matrix_MD_UA" in snakemake.input:
+        availability_MDUA = xr.open_dataarray(snakemake.input["availability_matrix_MD_UA"])
+        availability.loc[availability_MDUA.coords] = availability_MDUA
+        
     area = cutout.grid.to_crs(3035).area / 1e6
     area = xr.DataArray(area.values.reshape(cutout.shape),
                         [cutout.coords['y'], cutout.coords['x']])

scripts/determine_availability_matrix_MD_UA.py.ipynb

@@ -0,0 +1,126 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "1f7daec4",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import progressbar as pgb\n",
+    "import geopandas as gpd\n",
+    "import xarray as xr\n",
+    "import numpy as np\n",
+    "import functools\n",
+    "import atlite\n",
+    "import logging\n",
+    "import time\n",
+    "import matplotlib.pyplot as plt\n",
+    "from _helpers import mock_snakemake\n",
+    "\n",
+    "from _helpers import configure_logging\n",
+    "\n",
+    "logger = logging.getLogger(__name__)\n",
+    "\n",
+    "if __name__ == '__main__':\n",
+    "    if 'snakemake' not in globals():\n",
+    "        from _helpers import mock_snakemake\n",
+    "        snakemake = mock_snakemake('determine_availability_matrix_MD_UA', technology='solar')\n",
+    "    configure_logging(snakemake)\n",
+    "    pgb.streams.wrap_stderr()\n",
+    "\n",
+    "    nprocesses = snakemake.config['atlite'].get('nprocesses')\n",
+    "    noprogress = not snakemake.config['atlite'].get('show_progress', True)\n",
+    "    config = snakemake.config['renewable'][snakemake.wildcards.technology]\n",
+    "\n",
+    "    cutout = atlite.Cutout(snakemake.input['cutout'])\n",
+    "    regions = gpd.read_file(snakemake.input.regions).set_index('name').rename_axis('bus')\n",
+    "    buses = regions.index\n",
+    "\n",
+    "    excluder = atlite.ExclusionContainer(crs=3035, res=100)\n",
+    "\n",
+    "    corine = config.get(\"corine\", {})\n",
+    "    if \"grid_codes\" in corine:\n",
+    "        \n",
+    "        # Land cover codes to emulate CORINE results\n",
+    "        if snakemake.wildcards.technology == \"solar\":\n",
+    "            codes = [20, 30, 40, 50, 60, 90, 100]\n",
+    "        elif snakemake.wildcards.technology == \"onwind\":\n",
+    "            codes = [20, 30, 40, 60, 100, 111, 112, 113, 114, 115, 116, 121, 122, 123, 124, 125, 126]\n",
+    "        elif snakemake.wildcards.technology == \"offshore-ac\":\n",
+    "            codes = [80, 200]\n",
+    "        elif snakemake.wildcards.technology == \"offshore-dc\":\n",
+    "            codes = [80, 200]\n",
+    "        else:\n",
+    "            assert False, \"technology not supported\"\n",
+    "            \n",
+    "        excluder.add_raster(snakemake.input.copernicus, codes=codes, invert=True, crs=\"EPSG:4326\")\n",
+    "    if \"distance\" in corine and corine.get(\"distance\", 0.) > 0.:\n",
+    "        # Land cover codes to emulate CORINE results\n",
+    "        if snakemake.wildcards.technology == \"onwind\":\n",
+    "            codes = [50]\n",
+    "        else:\n",
+    "            assert False, \"technology not supported\"\n",
+    "            \n",
+    "        buffer = corine[\"distance\"]\n",
+    "        excluder.add_raster(snakemake.input.copernicus, codes=codes, buffer=buffer, crs=\"EPSG:4326\")\n",
+    "\n",
+    "    if \"max_depth\" in config:\n",
+    "        # lambda not supported for atlite + multiprocessing\n",
+    "        # use named function np.greater with partially frozen argument instead\n",
+    "        # and exclude areas where: -max_depth > grid cell depth\n",
+    "        func = functools.partial(np.greater,-config['max_depth'])\n",
+    "        excluder.add_raster(snakemake.input.gebco, codes=func, crs=4236, nodata=-1000)\n",
+    "\n",
+    "    if 'min_shore_distance' in config:\n",
+    "        buffer = config['min_shore_distance']\n",
+    "        excluder.add_geometry(snakemake.input.country_shapes, buffer=buffer)\n",
+    "\n",
+    "    if 'max_shore_distance' in config:\n",
+    "        buffer = config['max_shore_distance']\n",
+    "        excluder.add_geometry(snakemake.input.country_shapes, buffer=buffer, invert=True)\n",
+    "\n",
+    "    kwargs = dict(nprocesses=nprocesses, disable_progressbar=noprogress)\n",
+    "    if noprogress:\n",
+    "        logger.info('Calculate landuse availabilities...')\n",
+    "        start = time.time()\n",
+    "        availability = cutout.availabilitymatrix(regions, excluder, **kwargs)\n",
+    "        duration = time.time() - start\n",
+    "        logger.info(f'Completed availability calculation ({duration:2.2f}s)')\n",
+    "    else:\n",
+    "        availability = cutout.availabilitymatrix(regions, excluder, **kwargs)\n",
+    "\n",
+    "    # Limit results only to buses for UA and MD\n",
+    "    buses = regions.loc[regions[\"country\"].isin([\"UA\",\"MD\"])].index.values\n",
+    "    availability = availability.sel(bus=buses)\n",
+    "    \n",
+    "    # Save and plot for verification\n",
+    "    availability.to_netcdf(snakemake.output[\"availability_matrix\"])\n",
+    "    #availability.sum(dim=\"bus\").plot()\n",
+    "    #plt.title(technology)\n",
+    "    #plt.show()"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
+  },
+  "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.8.12"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}