pypsa-eur/notebooks/pypsa-eur.ipynb

2063 lines
64 KiB
Plaintext
Raw Normal View History

2017-10-11 23:51:24 +00:00
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Limitations\n",
"\n",
"*Please read and understand these carefully*\n",
"\n",
"While the benefit of an openly available, functional and partially\n",
"validated model of the European transmission system is high, many\n",
"approximations have been made due to missing data. In this section we\n",
"summarise the limitations of the dataset, both as a warning to the\n",
"user and as an encouragement to assist in improving the approximations.\n",
"\n",
"The grid data is based on a map of the ENTSO-E area\n",
"that is known to contain small distortions to improve\n",
"readability. Since the exact impedances of the lines are unknown,\n",
"approximations based on line lengths and standard line parameters\n",
"ignore specific conductoring choices for particular lines. There is no\n",
"openly available data on busbar configurations, switch locations,\n",
"transformers or reactive power compensation assets.\n",
"\n",
"Using Voronoi cells to aggregate load and generator data to\n",
"transmission network substations ignores the topology of the\n",
"underlying distribution network, meaning that assets may be connected\n",
"to the wrong substation. Assumptions have been made about the\n",
"distribution of load in each country proportional to population and\n",
"GDP that may not reflect local circumstances. Openly available data on\n",
"load time series may not correspond to the true vertical load and is not\n",
"spatially disaggregated; assuming, as we have done, that the load time\n",
"series shape is the same at each node within each country ignores local\n",
"differences.\n",
"\n",
"Wind, solar and small hydro, geothermal, marine and biomass power\n",
"plants are excluded from the dataset because of a lack of data\n",
"availability in many countries. Approximate distributions of wind and\n",
"solar plants in each country can be generated that are proportional to\n",
"the capacity factor at each location.\n",
"\n",
"The database of hydro-electric power plants does not include\n",
"plant-specific energy storage information, so that blanket values\n",
"based on country storage totals have been used. Inflow time series are\n",
"based on country-wide approximations, ignoring local topography and\n",
"basin drainage; in principle a full hydrological model should be used.\n",
"\n",
"Border connections and power flows to Russia, Belarus, Ukraine, Turkey\n",
"and Morocco have not been taken into account; islands which are not\n",
"connected to the main European system, such as Malta, Crete and\n",
"Cyprus, are also excluded from the model.\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Initialization"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import matplotlib as mpl\n",
"import matplotlib.pyplot as plt\n",
"%matplotlib inline\n",
"import seaborn as sns\n",
"sns.set_style('white')\n",
"\n",
"import pandas as pd\n",
"import numpy as np\n",
"import geopandas as gpd\n",
"\n",
"import scipy as sp, scipy.sparse\n",
"\n",
"import shapely.wkt\n",
"from shapely.geometry import Point\n",
"import shapely.prepared\n",
"\n",
"import os\n",
"from itertools import chain, count, ifilter\n",
"from six import itervalues, iteritems\n",
"from operator import attrgetter\n",
"\n",
"import xarray as xr"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import pypsa\n",
"\n",
"from vresutils import (shapes as vshapes, plot as vplot, landuse as vlanduse,\n",
" costdata2 as vcostdata2, hydro as vhydro)\n",
"from vresutils.graph import voronoi_partition_pts, BreadthFirstLevels\n",
"from vresutils.costdata import annuity, USD2013_to_EUR2013\n",
"\n",
"import powerplantmatching as pm\n",
"\n",
"import atlite"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Configuration"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"countries = ['AL', 'AT', 'BA', 'BE', 'BG', 'CH', 'CZ', 'DE', 'DK', 'EE', 'ES',\n",
" 'FI', 'FR', 'GB', 'GR', 'HR', 'HU', 'IE', 'IT', 'LT', 'LU', 'LV',\n",
" 'ME', 'MK', 'NL', 'NO', 'PL', 'PT', 'RO', 'RS', 'SE', 'SI', 'SK']"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Network topology"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Data files"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Buses"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"buses = (pd.read_csv(\"data/buses.csv\", quotechar=\"'\", true_values='t', false_values='f', dtype=dict(bus_id=\"str\"))\n",
" .set_index(\"bus_id\")\n",
" .rename(columns=dict(voltage='v_nom')))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"buses_pts = buses.pop('geometry').map(shapely.wkt.loads)\n",
"buses['x'] = buses_pts.map(attrgetter('x'))\n",
"buses['y'] = buses_pts.map(attrgetter('y'))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"buses['carrier'] = buses.pop('dc').map({True: 'DC', False: 'AC'})"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Transformers and Converters"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"transformers = (pd.read_csv(\"data/transformers.csv\", quotechar=\"'\", true_values='t', false_values='f',\n",
" dtype=dict(transformer_id='str', src_bus_id='str', dst_bus_id='str'))\n",
" .set_index('transformer_id')\n",
" .rename(columns=dict(src_bus_id='bus0', dst_bus_id='bus1', voltage='v_nom'))\n",
" .drop(['symbol'], axis=1))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# src_dc, dst_dc, src_voltage and dst_voltage contain redundant information,\n",
"# as can easily be verified\n",
"\n",
"assert ((transformers.src_dc == (transformers.bus0.map(buses.carrier) == 'DC'))\n",
" | transformers.src_dc.isnull()).all()\n",
"assert ((transformers.dst_dc == (transformers.bus1.map(buses.carrier) == 'DC'))\n",
" | transformers.dst_dc.isnull()).all()\n",
"assert ((transformers.src_voltage == transformers.bus0.map(buses.v_nom))\n",
" | transformers.src_voltage.isnull()).all()\n",
"assert ((transformers.dst_voltage == transformers.bus1.map(buses.v_nom))\n",
" | transformers.dst_voltage.isnull()).all()\n",
"\n",
"transformers.drop(['src_dc', 'dst_dc', 'src_voltage', 'dst_voltage'],\n",
" inplace=True, axis=1)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"There are transformers and converters in the dataset. Converters are those, which connect an AC to a DC bus."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"same_carrier = transformers.bus1.map(buses.carrier) == transformers.bus0.map(buses.carrier)\n",
"converters = transformers.loc[~ same_carrier]\n",
"transformers = transformers.loc[same_carrier]\n",
"del same_carrier"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### AC and DC lines"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"lines = (pd.read_csv(\"data/links.csv\", quotechar=\"'\", true_values='t', false_values='f',\n",
" dtype=dict(link_id='str', src_bus_id='str', dst_bus_id='str'))\n",
" .set_index('link_id')\n",
" .rename(columns=dict(src_bus_id='bus0', dst_bus_id='bus1', voltage='v_nom')))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"lines = lines.assign(length=lines['length_m'] / 1e3).drop('length_m', axis=1)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"dclines = lines.query('dc').drop('dc', axis=1)\n",
"aclines = lines.query('~ dc').drop('dc', axis=1)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Reduce dataset"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Remove duplicate lines"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"line_dups = pd.read_csv('data/line-dups.csv', dtype='str', skiprows=11)['line_id']\n",
"aclines.drop(line_dups, inplace=True)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"dclines = pd.DataFrame(dclines[~pd.DataFrame(dict(lat1=dclines.bus0.map(buses.y), lon1=dclines.bus0.map(buses.x),\n",
" lat2=dclines.bus1.map(buses.y), lon2=dclines.bus1.map(buses.x))).duplicated()])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Remove all but 220kV, 300kV and 380kV voltage levels"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"voltage_levels_to_remove = [110.0,\n",
" 150.0,\n",
" 132.0,\n",
" 750.0]"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"buses_with_v_nom_to_keep_b = ~ buses.v_nom.isin(voltage_levels_to_remove)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"transformers_with_v_nom_to_keep_b = transformers.bus0.map(buses_with_v_nom_to_keep_b) & transformers.bus1.map(buses_with_v_nom_to_keep_b)\n",
"converters_with_v_nom_to_keep_b = converters.bus0.map(buses_with_v_nom_to_keep_b) & converters.bus1.map(buses_with_v_nom_to_keep_b)\n",
"aclines_with_v_nom_to_keep_b = aclines.bus0.map(buses_with_v_nom_to_keep_b) & aclines.bus1.map(buses_with_v_nom_to_keep_b)\n",
"dclines_with_v_nom_to_keep_b = dclines.bus0.map(buses_with_v_nom_to_keep_b) & dclines.bus1.map(buses_with_v_nom_to_keep_b)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Remove everything not in any of the configured countries"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"europe_shape = vshapes.country_cover(countries)\n",
"europe_shape_exterior = shapely.geometry.Polygon(shell=europe_shape.exterior) # no holes"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"vplot.shapes(europe_shape_exterior, colour='blue').set_zorder(-1)\n",
"vplot.draw_basemap()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import shapely.prepared\n",
"europe_shape_exterior_prepped = shapely.prepared.prep(europe_shape_exterior)\n",
"buses_in_europe_b = buses_pts.map(lambda p: europe_shape_exterior_prepped.contains(p))\n",
"del europe_shape_exterior_prepped\n",
"#buses = buses.loc[buses_in_europe]"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"transformers_in_europe_b = transformers.bus0.map(buses_in_europe_b) | transformers.bus1.map(buses_in_europe_b)\n",
"converters_in_europe_b = converters.bus0.map(buses_in_europe_b) | converters.bus1.map(buses_in_europe_b)\n",
"aclines_in_europe_b = aclines.bus0.map(buses_in_europe_b) | aclines.bus1.map(buses_in_europe_b)\n",
"dclines_in_europe_b = dclines.bus0.map(buses_in_europe_b) | dclines.bus1.map(buses_in_europe_b)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from functools import reduce\n",
"\n",
"buses_in_europe_b |= pd.Series(\n",
" buses.index.isin(reduce(set.union,\n",
" (x[i]\n",
" for x in (converters.loc[converters_in_europe_b],\n",
" transformers.loc[transformers_in_europe_b],\n",
" aclines.loc[aclines_in_europe_b],\n",
" dclines.loc[dclines_in_europe_b])\n",
" for i in ('bus0', 'bus1')),\n",
" set())),\n",
" index=buses.index)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Perform prepared removals"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"buses = pd.DataFrame(buses.loc[buses_in_europe_b & buses_with_v_nom_to_keep_b])\n",
"converters = pd.DataFrame(converters.loc[converters_in_europe_b & converters_with_v_nom_to_keep_b])\n",
"transformers = pd.DataFrame(transformers.loc[transformers_in_europe_b & transformers_with_v_nom_to_keep_b])\n",
"aclines = pd.DataFrame(aclines.loc[aclines_in_europe_b & aclines_with_v_nom_to_keep_b])\n",
"dclines = pd.DataFrame(dclines.loc[dclines_in_europe_b & dclines_with_v_nom_to_keep_b])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Adjust parameters"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def find_closest_bus(pos):\n",
" if (not hasattr(find_closest_bus, 'kdtree')) or len(find_closest_bus.kdtree.data) != len(buses.index):\n",
" find_closest_bus.kdtree = sp.spatial.cKDTree(buses.loc[:,[\"x\", \"y\"]].values)\n",
" return buses.index[find_closest_bus.kdtree.query(pos)[1]]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Transformers"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"transformers = transformers.append(\n",
" pd.DataFrame([(str(transformers.index.astype(int).max()+1),\n",
" find_closest_bus((22.6, 53.9)),\n",
" find_closest_bus((22.5, 53.9)))],\n",
" columns=['index', 'bus0', 'bus1']).set_index('index')\n",
")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"transformers[\"x\"] = 0.1\n",
"transformers[\"s_nom\"] = 2000\n",
"transformers['type'] = \"\""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Converters + DC-Lines ->Links"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"converters = converters.append(\n",
" pd.DataFrame([\n",
" (i, find_closest_bus(pos1), find_closest_bus(pos2))\n",
" for i, (pos1, pos2) in enumerate([((6.8, 59.6), (6.65, 59.5)),\n",
" ((16.7, 57.6), (16.7, 57.4)),\n",
" ((1.3, 51.2), (1.1, 51.3))],\n",
" start=converters.index.astype(int).max()+1)\n",
" ], columns=['index', 'bus0', 'bus1']).set_index('index')\n",
")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"converters = converters.groupby('bus0', as_index=False).apply(\n",
" lambda x: x.iloc[0] if len(x) == 1 else x.loc[x['bus1'].map(buses['v_nom']).idxmax()]\n",
")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"converters_ac_to_dc_b = (converters.bus0.map(buses['carrier']) == 'AC') & (converters.bus1.map(buses['carrier']) == 'DC')\n",
"dc_to_ac_map = (pd.Series(converters.loc[converters_ac_to_dc_b, 'bus0'].values, index=converters.loc[converters_ac_to_dc_b, 'bus1'].values)\n",
" .append(pd.Series(converters.loc[~converters_ac_to_dc_b, 'bus1'].values, index=converters.loc[~converters_ac_to_dc_b, 'bus0'].values)))\n",
"missing_dc_buses_i = buses[buses['carrier'] == 'DC'].index.difference(dc_to_ac_map.index)\n",
"dc_to_ac_map = dc_to_ac_map.append(pd.Series(missing_dc_buses_i, index=missing_dc_buses_i))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Drop unneeded buses"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"buses.drop(dc_to_ac_map.index[dc_to_ac_map.index.values != dc_to_ac_map.values], inplace=True)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Rewire dclines"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"dclines['bus0'] = dclines['bus0'].map(dc_to_ac_map)\n",
"dclines['bus1'] = dclines['bus1'].map(dc_to_ac_map)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"dclines = dclines.append(\n",
" pd.DataFrame([(str(dclines.index.astype(int).max()+1),\n",
" 500,\n",
" find_closest_bus((23.9, 54.4)),\n",
" find_closest_bus((24., 54.4)),\n",
" vshapes.haversine([(23.9, 54.4), (24., 54.4)]))],\n",
" columns=['index', 'p_nom', 'bus0', 'bus1', 'length']).set_index('index')\n",
")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"dclines.drop('v_nom', axis=1, inplace=True)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### DC lines"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"DC line capacities from [List of HVDC Projects on Wikipedia](https://en.wikipedia.org/wiki/List_of_HVDC_projects)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"dclines['p_nom'] = pd.read_csv('data/dclines_p-nom.csv', names=['name', 'p_nom'], dtype={'name': str, 'p_nom': float}, skiprows=1).set_index('name')\n",
"dclines['p_min_pu'] = -1."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### AC-Lines"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Split AC-Lines with multiple voltages"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def namer(string, start=0): return (string.format(x) for x in count(start=start))\n",
"busname = namer(\"M{:02}\")\n",
"trafoname = namer(\"M{:02}\")\n",
"linename = namer(\"M{:02}\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def find_or_add_lower_v_nom_bus(bus, v_nom):\n",
" candidates = transformers.loc[(transformers.bus1 == bus) & (transformers.bus0.map(buses.v_nom) == v_nom), 'bus0']\n",
" if len(candidates):\n",
" return candidates.iloc[0]\n",
" new_bus = next(busname)\n",
" buses.loc[new_bus] = pd.Series({'v_nom': v_nom, 'symbol': 'joint', 'carrier': 'AC',\n",
" 'x': buses.at[bus, 'x'], 'y': buses.at[bus, 'y'],\n",
" 'under_construction': buses.at[bus, 'under_construction']})\n",
" transformers.loc[next(trafoname)] = pd.Series({'bus0': new_bus, 'bus1': bus, 'x': 0.1, 's_nom': 2000, 'type': \"\"})\n",
" return new_bus"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"voltage_levels = aclines.v_nom.unique()\n",
"\n",
"for line in aclines.tags.str.extract(r'\"text_\"=>\"\\(?(\\d+)\\+(\\d+).*?\"', expand=True).dropna().itertuples():\n",
" v_nom = int(line._2)\n",
" if aclines.at[line.Index, 'circuits'] > 1:\n",
" aclines.at[line.Index, 'circuits'] -= 1\n",
"\n",
" if v_nom in voltage_levels:\n",
" bus0 = find_or_add_lower_v_nom_bus(aclines.at[line.Index, 'bus0'], v_nom)\n",
" bus1 = find_or_add_lower_v_nom_bus(aclines.at[line.Index, 'bus1'], v_nom)\n",
" aclines.loc[next(linename)] = pd.Series(dict(chain(iteritems({'bus0': bus0, 'bus1': bus1, 'v_nom': v_nom, 'circuits': 1}),\n",
" iteritems({k: aclines.at[line.Index, k]\n",
" for k in {'underground', 'under_construction',\n",
" 'tags', 'geometry', 'length'}}))))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Electrical parameters"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Apply standard line-types from https://pypsa.org/doc/components.html#line-types"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# for v_nom, linetype in [(220., \"Al/St 240/40 2-bundle 220.0\"),\n",
"# (300., \"Al/St 240/40 3-bundle 300.0\"),\n",
"# (380., \"Al/St 240/40 4-bundle 380.0\")]:\n",
"# aclines.loc[aclines[\"v_nom\"] == v_nom, 'type'] = linetype"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from collections import namedtuple\n",
"\n",
"LineParam = namedtuple(\"LineParam\", [\"v_nom\", \"wires\", \"r\", \"x\", \"c\", \"i\"])\n",
"for p in (LineParam(v_nom=220., wires=2.0, r=0.06, x=0.301, c=12.5, i=1.29),\n",
" LineParam(v_nom=380., wires=4.0, r=0.03, x=0.246, c=13.8, i=2.58),\n",
" LineParam(v_nom=300., wires=3.0, r=0.04, x=0.265, c=13.2, i=1.935)):\n",
" ls = aclines[\"v_nom\"] == p.v_nom\n",
" length = aclines.loc[ls, \"length\"]\n",
" circuits = aclines.loc[ls, \"circuits\"]\n",
" aclines.loc[ls, \"r\"] = length * p.r / circuits\n",
" aclines.loc[ls, \"x\"] = length * p.x / circuits\n",
" aclines.loc[ls, \"b\"] = length * (2*np.pi*50*1e-9) * p.c * circuits\n",
" aclines.loc[ls, \"s_nom\"] = np.sqrt(3.) * p.v_nom * p.i * circuits"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"aclines['type'] = \"\""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## PyPSA model"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"network = pypsa.Network()\n",
"network.name = \"PyPSA-Europe\"\n",
"network"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"pypsa.io.import_components_from_dataframe(network, buses, \"Bus\")\n",
"pypsa.io.import_components_from_dataframe(network, transformers, \"Transformer\")\n",
"pypsa.io.import_components_from_dataframe(network, dclines, 'Link')\n",
"pypsa.io.import_components_from_dataframe(network, aclines, 'Line')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Connect close buses (<= 1km distance)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def connect_close_buses(network, radius=1.):\n",
" adj = network.graph([\"Line\", \"Transformer\", \"Link\"]).adj\n",
" \n",
" n_lines_added = 0\n",
" n_transformers_added = 0\n",
" ac_buses = network.buses[network.buses.carrier == 'AC']\n",
"\n",
" for i,u in enumerate(ac_buses.index):\n",
"\n",
" vs = ac_buses[[\"x\",\"y\"]].iloc[i+1:]\n",
" distance_km = pypsa.geo.haversine(vs, ac_buses.loc[u,[\"x\",\"y\"]])\n",
" \n",
" for j,v in enumerate(vs.index):\n",
" km = distance_km[j,0]\n",
" \n",
" if km < radius:\n",
" if u in adj[v]:\n",
" continue\n",
" #print(u,v,ac_buses.at[u,\"v_nom\"], ac_buses.at[v,\"v_nom\"],km)\n",
" \n",
" if ac_buses.at[u,\"v_nom\"] != ac_buses.at[v,\"v_nom\"]:\n",
" network.add(\"Transformer\",\"extra_trafo_{}_{}\".format(u,v),s_nom=2000,bus0=u,bus1=v,x=0.1)\n",
" n_transformers_added += 1\n",
" else:\n",
" network.add(\"Line\",\"extra_line_{}_{}\".format(u,v),s_nom=4000,bus0=u,bus1=v,x=0.1)\n",
" n_lines_added += 1\n",
" \n",
" print(\"Added {} lines and {} transformers to connect buses less than {} km apart.\"\n",
" .format(n_lines_added, n_transformers_added, radius))\n",
"\n",
"connect_close_buses(network)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Remove all connected components with less than 10 buses"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"network.determine_network_topology()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"sub_network_sizes = network.buses.groupby('sub_network').size()\n",
"subs_to_remove = sub_network_sizes.index[sub_network_sizes < 10]"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"print(\"Removing {} small sub_networks (synchronous zones) with less than 10 buses. In total {} buses.\"\n",
" .format(len(subs_to_remove), network.buses.sub_network.isin(subs_to_remove).sum()))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"network = network[~network.buses.sub_network.isin(subs_to_remove)]"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"fig, ax = plt.subplots()\n",
"voltage_colors = {132.0: 'blue', 220.0: 'green', 300.0: 'orange', 380.0: 'red'}\n",
"network.plot(bus_sizes=5,\n",
" bus_colors=network.buses['v_nom'].map(voltage_colors),\n",
" line_colors=pd.concat(dict(Line=network.lines['v_nom'].map(voltage_colors),\n",
" Link=pd.Series('violet', index=dclines.index))))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Export network topology"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"network.export_to_csv_folder('results/01-network-topology')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Bus regions"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"country_shapes = vshapes.countries(countries, minarea=0.1, tolerance=0.01, add_KV_to_RS=True)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"offshore_shapes = vshapes.eez(countries, tolerance=0.01)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"On- and offshore area should be distributed to all buses"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"vplot.shapes(offshore_shapes)\n",
"network.plot(bus_sizes=20)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def buses_in_shape(shape):\n",
" shape = shapely.prepared.prep(shape)\n",
" return pd.Series(\n",
" np.fromiter((shape.contains(Point(x, y))\n",
" for x, y in network.buses.loc[:,[\"x\", \"y\"]].values),\n",
" dtype=bool, count=len(network.buses)),\n",
" index=network.buses.index\n",
" )"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"As there are multiple voltage levels and AC and DC buses, we preferable want to attach load and onshore generation to low voltage AC buses, while the offshore generation should go to high voltage AC buses, if there is no offshore bus built yet."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"network.buses['substation'] = network.buses.symbol.str.contains('substation', case=False)\n",
"\n",
"def prefer_ac_voltage(x, which):\n",
" index = x.index\n",
" if len(index) == 1:\n",
" return pd.Series(index, index)\n",
" ac = x['carrier'] == 'AC'\n",
" if ac.sum() > 0:\n",
" x = x.loc[ac]\n",
" key = (x.index[0]\n",
" if x['v_nom'].isnull().all()\n",
" else getattr(x['v_nom'], 'idx' + which)())\n",
" return pd.Series(key, index)\n",
"\n",
"gb = network.buses.loc[lambda df: df.substation].groupby(['x', 'y'], as_index=False, group_keys=False, sort=False)\n",
"bus_map_low = gb.apply(prefer_ac_voltage, 'min') #.reindex(network.buses.index)\n",
"buses_onshore_has_region_b = (bus_map_low == bus_map_low.index).reindex(network.buses.index, fill_value=False)\n",
"bus_map_high = gb.apply(prefer_ac_voltage, 'max')\n",
"buses_onshore_for_offshore_has_region_b = (bus_map_high == bus_map_high.index).reindex(network.buses.index, fill_value=False)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Split countries into voronoi cells and attach them to the buses as region_onshore and region_offshore entries"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"%%time\n",
"network.buses.drop(['region_onshore', 'region_offshore'], inplace=True, axis=1, errors='ignore')\n",
"\n",
"for country in countries:\n",
" onshore_shape = country_shapes[country]\n",
"\n",
" # onshore\n",
" buses_onshore_b = buses_in_shape(onshore_shape)\n",
" network.buses.loc[buses_onshore_b, 'country'] = country \n",
"\n",
" buses_onshore_locations = network.buses.loc[buses_onshore_b & buses_onshore_has_region_b, [\"x\", \"y\"]]\n",
" network.buses.loc[buses_onshore_locations.index, 'region_onshore'] = voronoi_partition_pts(buses_onshore_locations.values, onshore_shape)\n",
"\n",
" # offshore\n",
" if country not in offshore_shapes: continue\n",
" offshore_shape = offshore_shapes[country]\n",
"\n",
" buses_offshore_b = buses_in_shape(offshore_shape)\n",
" network.buses.loc[buses_offshore_b, 'country'] = country\n",
" \n",
" buses_offshore_locations = network.buses.loc[buses_offshore_b | (buses_onshore_b & buses_onshore_for_offshore_has_region_b), [\"x\", \"y\"]]\n",
" buses_offshore_regions = pd.Series(voronoi_partition_pts(buses_offshore_locations.values, offshore_shape),\n",
" index=buses_offshore_locations.index)\n",
" buses_offshore_regions = buses_offshore_regions.loc[buses_offshore_regions.map(attrgetter('area')) > 1e-2]\n",
" network.buses.loc[buses_offshore_regions.index, 'region_offshore'] = buses_offshore_regions"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"onshorebuses_i = network.buses.index[network.buses['region_onshore'].notnull()]\n",
"offshorebuses_i = network.buses.index[network.buses['region_offshore'].notnull()]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Some buses don't lie in any country. Just add them to the nearest (hop-wise) country, but don't give them any load or generation, i.e. they will not be part of onshorebuses or offshorebuses."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"buses_nan_b = network.buses.country.isnull()\n",
"graph = network.graph()\n",
"network.buses.loc[buses_nan_b, 'country'] = \\\n",
" [(next(ifilter(len, map(lambda x: network.buses.loc[x, 'country'].dropna(),\n",
" BreadthFirstLevels(graph, [b]))))\n",
" .value_counts().index[0])\n",
" for b in network.buses.index[buses_nan_b]]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Compute area of each region"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"network.buses.drop(['area', 'area_offshore'], axis=1, inplace=True, errors='ignore')\n",
"\n",
"on_b = network.buses['region_onshore'].notnull()\n",
"network.buses.loc[on_b, 'area'] = network.buses.loc[on_b, 'region_onshore'].map(vshapes.area)\n",
"\n",
"off_b = network.buses['region_offshore'].notnull()\n",
"network.buses.loc[off_b, 'area_offshore'] = network.buses.loc[off_b, 'region_offshore'].map(vshapes.area)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Illustrate regions of DE"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"busesDE = network.buses.query('country == \"DE\"')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Offshore"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"plt.scatter(*busesDE.loc[busesDE.region_offshore.notnull(), ['x', 'y']].values.T).set_zorder(2)\n",
"vplot.shapes(busesDE.loc[busesDE.region_offshore.notnull(), 'region_offshore'])\n",
"vplot.shapes([offshore_shapes['DE']], outline=True, colour='r')\n",
"vplot.draw_basemap()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Onshore"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"plt.scatter(*busesDE.loc[busesDE.region_onshore.notnull(), ['x', 'y']].values.T).set_zorder(2)\n",
"vplot.shapes(busesDE.loc[busesDE.region_onshore.notnull(), 'region_onshore'])\n",
"vplot.draw_basemap()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Generation"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"network.set_snapshots(pd.date_range(start='2011', freq='h', periods=5*8760 + 24))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Costs"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"costs = vcostdata2.get_full_cost_CO2('diw2030', 0.)\n",
"costs"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Carriers"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"network.carriers.drop(network.carriers.index, inplace=True)\n",
"\n",
"carriers = pd.DataFrame([('onwind', 'WON'),\n",
" ('offwind', 'WOFF'),\n",
" ('solar', 'S'),\n",
" ('PHS', 'Hy'),\n",
" ('hydro', 'Hy'),\n",
" ('ror', 'Hy'),\n",
" #('H2', 'H'),\n",
" #('battery', 'B'),\n",
" ('OCGT', 'GO'),\n",
" ('CCGT', 'GC'),\n",
" ('coal', 'C'),\n",
" ('nuclear', 'N'),\n",
" ('lignite', 'L')],\n",
" columns=['name', 'short_name']).set_index('name')\n",
"carriers['co2_emissions'] = costs.loc[carriers.index, 'CO2intensity'].fillna(0.)\n",
"\n",
"network.import_components_from_dataframe(carriers, \"Carrier\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Conventional Generation "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"powerplants = pm.MATCHED_dataset(subsume_uncommon_fueltypes=True, include_unavailables=True)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Fuel types"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Nearly all of the natural gas systems are CCGT"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"conventionals = pd.DataFrame(powerplants.loc[(~powerplants.Fueltype.isin(['Wind', 'Hydro', 'Solar']))\n",
" & powerplants.lon.notnull() & powerplants.lat.notnull()])"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"conventionals.loc[conventionals.Fueltype == 'Natural Gas', 'Technology'].value_counts(dropna=False)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"conventionals.loc[(conventionals.Fueltype == 'Natural Gas') & conventionals.Technology.str.contains('OCGT|Open'), 'Fueltype'] = 'OCGT'\n",
"conventionals.loc[(conventionals.Fueltype == 'Natural Gas') & (~conventionals.Technology.str.contains('OCGT|Open').fillna(False)), 'Fueltype'] = 'CCGT'"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"conventionals = conventionals.loc[conventionals.Fueltype.isin(['Hard Coal', 'Lignite', 'Nuclear', 'OCGT', 'CCGT', 'Hydro'])]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Add costs, efficiency and a unique name"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"conventionals.Fueltype.replace({'Hard Coal': 'coal', 'Nuclear': 'nuclear', 'Lignite': 'lignite'}, inplace=True)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"for ft in conventionals.Fueltype.unique():\n",
" conv_ft_b = conventionals.Fueltype == ft\n",
" conventionals.loc[conv_ft_b, 'marginal_cost'] = costs.at[ft, 'marginal']\n",
" conventionals.loc[conv_ft_b, 'capital_cost'] = costs.at[ft, 'capital']\n",
" conventionals.loc[conv_ft_b, 'efficiency'] = costs.at[ft, 'efficiency']\n",
" conventionals.loc[conv_ft_b, 'name'] = carriers.at[ft, 'short_name'] + pd.RangeIndex(conv_ft_b.sum()).astype(str)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Locate buses"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"kdtree = sp.spatial.cKDTree(network.buses.loc[onshorebuses_i, ['x','y']].values)\n",
"conventionals.loc[:, 'bus'] = onshorebuses_i[kdtree.query(conventionals[['lon','lat']].values)[1]]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Import into pypsa"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"conventionals = conventionals.rename(columns={'Capacity': 'p_nom', 'Fueltype': 'carrier'})\n",
"conventionals = conventionals[['bus', 'name', 'p_nom', 'carrier', 'marginal_cost', 'capital_cost', 'efficiency']].set_index('name')"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"network.import_components_from_dataframe(conventionals, 'Generator')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Hydro Generation"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"hydro = pd.DataFrame(powerplants.loc[lambda df: (df['Fueltype'] == 'Hydro') & df['lon'].notnull() & df['lat'].notnull()])\n",
"hydro.set_index(carriers.at['hydro', 'short_name'] + pd.RangeIndex(len(hydro)).astype(str), drop=False, inplace=True)\n",
"\n",
"hydro.loc[:, 'bus'] = onshorebuses_i[kdtree.query(hydro[['lon', 'lat']].values)[1]]"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"hydro = hydro.rename(columns={'Capacity':'p_nom', 'Technology':'technology'})\n",
"hydro = hydro.loc[hydro.technology.notnull(), ['bus', 'p_nom', 'technology']]"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"hydro_capas = vhydro.get_hydro_capas()\n",
"hydro_capas.loc[hydro_capas['E_store[TWh]'] < 0.2, 'E_store[TWh]'] = 0.2"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Inflow"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def normed(x):\n",
" return x/x.sum()\n",
"\n",
"country_shapes = pd.Series(country_shapes)\n",
"country_shapes.index.rename('countries', inplace=True)\n",
"country_shapes = country_shapes.reindex(countries)\n",
"\n",
"cutout_2011_2016 = atlite.Cutout('europe-2011-2016')\n",
"inflow = cutout_2011_2016.runoff(shapes=country_shapes,\n",
" smooth=True,\n",
" lower_threshold_quantile=True,\n",
" normalize_using_yearly=vhydro.get_eia_annual_hydro_generation().reindex(columns=countries))\n",
"\n",
"hydro_w_inflow = hydro.query(\"technology != 'Pumped Storage'\")\n",
"hydro_cntry = hydro_w_inflow.bus.map(network.buses.country)\n",
"\n",
"inflow_t = \\\n",
"(inflow.sel(countries=hydro_cntry.values)\n",
" .rename({'countries': 'name'})\n",
" .assign_coords(name=hydro_w_inflow.index)\n",
" .transpose('time', 'name')\n",
" .to_pandas()\n",
" .multiply(hydro_w_inflow.p_nom.groupby(hydro_cntry).transform(normed), axis=1)\n",
" .divide(hydro_w_inflow.p_nom, axis=1))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### RoR"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"hydro_ror = hydro.loc[hydro.technology == 'Run-Of-River']\n",
"\n",
"hydro_ror = hydro_ror.assign(\n",
" carrier=\"ror\",\n",
" efficiency=costs.at['ror', 'efficiency'],\n",
" weight=hydro_ror['p_nom']\n",
")\n",
"\n",
"pypsa.io.import_components_from_dataframe(network, hydro_ror, 'Generator')\n",
"\n",
"rorpotential = inflow_t.loc[:,hydro_ror.index].where(lambda df: df<=1., other=1.)\n",
"pypsa.io.import_series_from_dataframe(network, rorpotential, 'Generator', 'p_max_pu')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### PHS"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"hydro_phs = hydro.loc[hydro.technology == 'Pumped Storage']\n",
"\n",
"hydro_phs = hydro_phs.assign(\n",
" carrier='PHS',\n",
" max_hours=6,\n",
" p_max_pu=1.,\n",
" p_min_pu=-1.,\n",
" efficiency_store=np.sqrt(costs.at['PHS','efficiency']),\n",
" efficiency_dispatch=np.sqrt(costs.at['PHS','efficiency']),\n",
" cyclic_state_of_charge=True\n",
")\n",
"\n",
"network.import_components_from_dataframe(hydro_phs, 'StorageUnit')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Reservoir"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"hydro_reservoir = hydro.loc[hydro.technology == 'Reservoir']\n",
"\n",
"hydro_reservoir_max_hours = hydro_capas['E_store[TWh]']*1e6/hydro_reservoir.p_nom.groupby(hydro_cntry).sum()\n",
"\n",
"hydro_reservoir = hydro_reservoir.assign(\n",
" carrier=\"hydro\",\n",
" max_hours=hydro_cntry.loc[hydro_reservoir.index].map(hydro_reservoir_max_hours),\n",
" p_max_pu=1., #dispatch\n",
" p_min_pu=0.,\n",
" efficiency_dispatch=costs.at['hydro','efficiency'],\n",
" efficiency_store=0,\n",
" cyclic_state_of_charge=True,\n",
")\n",
"\n",
"network.import_components_from_dataframe(hydro_reservoir, 'StorageUnit')\n",
"\n",
"reservoirpotential = inflow_t.loc[:,hydro_reservoir.index]\n",
"pypsa.io.import_series_from_dataframe(network, reservoirpotential, 'StorageUnit', 'inflow')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Renewables"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"cutout = atlite.Cutout('europe-2011-2016')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Landuse potentials"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"%%time\n",
"onshorepotentials = xr.DataArray(vlanduse.windonshorepotentials(cutout, cushion_factor=0.3, distance=1000),\n",
" [cutout.coords['y'], cutout.coords['x']])\n",
"offshorepotentials = xr.DataArray(vlanduse.windoffshorepotentials(cutout, cushion_factor=0.3),\n",
" [cutout.coords['y'], cutout.coords['x']])\n",
"solarpotentials = xr.DataArray(vlanduse.solarpotentials(cutout),\n",
" [cutout.coords['y'], cutout.coords['x']])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Indicatormatrices and capacities"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"onshoreindicatormatrix = cutout.indicatormatrix(network.buses.loc[onshorebuses_i, 'region_onshore'])"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"offshoreindicatormatrix = cutout.indicatormatrix(network.buses.loc[offshorebuses_i, 'region_offshore'])"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Remove raster cells from indicatormatrix with very low potential (incl. those beyond 50m depth)\n",
"ocean_depth_cutoff = 50\n",
"offshorepotentials.values[(cutout.meta['height'] < - ocean_depth_cutoff).transpose(*offshorepotentials.dims)] = 0.\n",
"\n",
"offshoreindicatormatrix[:,(offshorepotentials < 1e-4).stack(spatial=('y', 'x')).values] = 0.\n",
"\n",
"keep_offshore_b = np.asarray(offshoreindicatormatrix.sum(axis=1)).squeeze() > 0.\n",
"network.buses.loc[offshorebuses_i[~keep_offshore_b], 'region_offshore'] = np.nan\n",
"network.buses.loc[offshorebuses_i[~keep_offshore_b], 'area_offshore'] = np.nan\n",
"offshoreindicatormatrix = offshoreindicatormatrix[keep_offshore_b]\n",
"offshorebuses_i = offshorebuses_i[keep_offshore_b]"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"onshore_i = ('WON' + pd.RangeIndex(len(onshorebuses_i)).astype(str)).rename(\"name\")\n",
"offshore_i = ('WOF' + pd.RangeIndex(len(offshorebuses_i)).astype(str)).rename(\"name\")\n",
"solar_i = ('S' + pd.RangeIndex(len(onshorebuses_i)).astype(str)).rename(\"name\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Generation time-series"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"onshorewindturbine ='Vestas_V112_3MW'\n",
"offshorewindturbine='NREL_ReferenceTurbine_5MW_offshore'\n",
"solarpanel = 'KANENA'\n",
"orientation='latitude_optimal'"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"%%time\n",
"onshorewindcapacityfactor = cutout.wind(turbine=onshorewindturbine, smooth=True, capacity_factor=True)\n",
"onshorewindlayout = onshorewindcapacityfactor * onshorepotentials\n",
"onshore, onshore_weights = cutout.wind(turbine=onshorewindturbine, smooth=True,\n",
" matrix=onshoreindicatormatrix, index=onshore_i,\n",
" layout=onshorewindlayout,\n",
" per_unit=True, return_units=True)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"%%time\n",
"offshorewindcapacityfactor = cutout.wind(turbine=offshorewindturbine, smooth=True, capacity_factor=True)\n",
"offshorewindlayout = offshorewindcapacityfactor * offshorepotentials\n",
"offshore, offshore_weights = cutout.wind(turbine=offshorewindturbine, smooth=True,\n",
" matrix=offshoreindicatormatrix, index=offshore_i,\n",
" layout=offshorewindlayout,\n",
" per_unit=True, return_units=True)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"%%time\n",
"solarlayout = cutout.pv(panel=solarpanel, orientation=orientation, capacity_factor=True) * solarpotentials\n",
"solar, solar_weights = cutout.pv(panel=solarpanel, orientation=orientation,\n",
" matrix=onshoreindicatormatrix, index=solar_i,\n",
" layout=solarlayout,\n",
" per_unit=True, return_units=True)\n",
"\n",
"# solar capacity factors in REatlas are too high compared to other\n",
"# studies; we correct this by applying a constant 'inverter\n",
"# inefficiency' to the p_max_pu timeseries; comparing with\n",
"# Pietzker+(2014) http://dx.doi.org/10.1016/j.apenergy.2014.08.011\n",
"# the produced-power weighted average of capacity factor ratios is 1.2683\n",
"solar /= 1.2683"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Determine maximal capacity per Voronoi cell"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def p_nom_max(potential, layout, indicatormatrix, weights):\n",
" weights = weights.to_series()\n",
" relativepotential = (potential / layout).stack(spatial=('y', 'x')).values\n",
" return pd.Series([np.nanmin(relativepotential[row.nonzero()[1]])\n",
" for row in indicatormatrix.tocsr()],\n",
" index=weights.index) * weights"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"onshore_p_nom_max = p_nom_max(onshorepotentials, onshorewindlayout, onshoreindicatormatrix, onshore_weights)\n",
"offshore_p_nom_max = p_nom_max(offshorepotentials, offshorewindlayout, offshoreindicatormatrix, offshore_weights)\n",
"solar_p_nom_max = p_nom_max(solarpotentials, solarlayout, onshoreindicatormatrix, solar_weights)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Illustrate German capacity factors and available capacity"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"onshore_de = pd.Series(onshorebuses_i, onshore_i).map(network.buses.country) == 'DE'\n",
"offshore_de = pd.Series(offshorebuses_i, offshore_i).map(network.buses.country) == 'DE'\n",
"solar_de = pd.Series(onshorebuses_i, solar_i).map(network.buses.country) == 'DE'\n",
"bins = np.r_[0.05:0.5:0.01]\n",
"#bins = np.r_[0.08:0.15:0.005]\n",
"plt.figure()\n",
"pd.concat(\n",
" dict(onshore=onshore_p_nom_max.loc[onshore_de].groupby(pd.cut(onshore.mean('time').to_series().loc[onshore_de], bins)).sum(),\n",
" offshore=offshore_p_nom_max.loc[offshore_de].groupby(pd.cut(offshore.mean('time').to_series().loc[offshore_de], bins)).sum(),\n",
" solar=solar_p_nom_max.loc[solar_de].groupby(pd.cut(solar.mean('time').to_series().loc[solar_de], bins)).sum()),\n",
" axis=1\n",
").fillna(0.).pipe(lambda df: df/1e3).plot(marker='o', rot=90)\n",
"plt.ylabel(\"Max Capacity / GW\")\n",
"plt.tight_layout()\n",
"\n",
"if False:\n",
" plt.savefig('capacityfactors_de.pdf')\n",
"#plt.savefig('capacityfactors_solar_de.pdf')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Plot onshore and offshore potentials and voronoi cells"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"grid_cells = pd.Series(cutout.grid_cells())\n",
"cell_areas = vlanduse._cutout_cell_areas(cutout)\n",
"onshorede = np.where(network.buses.loc[onshorebuses_i, 'country'] == 'DE')\n",
"offshorede = np.where(network.buses.loc[offshorebuses_i, 'country'] == 'DE')"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"onlayout = pd.Series((onshorewindlayout.transpose('y', 'x').values / cell_areas * 8760 / 1e3).ravel())\n",
"offlayout = pd.Series((offshorewindlayout.transpose('y', 'x').values / cell_areas * 8760 / 1e3).ravel())\n",
"solarlayout = pd.Series((solarlayout.transpose('y', 'x').values / cell_areas * 8760 / 1e3).ravel())"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# onbuses_i = np.concatenate([onshoreindicatormatrix[busno].nonzero()[1]\n",
"# for busno in busnos]\n",
"windpotential_per_grid_cell = (\n",
" onlayout * pd.Series(np.asarray(onshoreindicatormatrix[onshorede].sum(axis=0)).squeeze())\n",
" + offlayout * pd.Series(np.asarray(offshoreindicatormatrix[offshorede].sum(axis=0)).squeeze())\n",
").loc[lambda s: s>0]"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"solarpotential_per_grid_cell = (\n",
" solarlayout * pd.Series(np.asarray(onshoreindicatormatrix[onshorede].sum(axis=0)).squeeze())\n",
").loc[lambda s: s>0]"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"fig, ax = plt.subplots(figsize=(10, 6))\n",
"\n",
"ax.set_title('DE Wind average power density [GWh/a/km$^2$]')\n",
"\n",
"sh = vplot.shapes(grid_cells, windpotential_per_grid_cell, cmap='Blues',\n",
" norm=mpl.colors.Normalize(vmax=10), ax=ax)\n",
"vplot.shapes(network.buses.loc[onshorebuses_i[onshorede], 'region_onshore'],\n",
" outline=True, colour='black', linestyle='-', ax=ax)\n",
"vplot.shapes(network.buses.loc[offshorebuses_i[offshorede], 'region_offshore'],\n",
" outline=True, colour='black', linestyle='-', ax=ax)\n",
"fig.colorbar(sh, extend='max')\n",
"\n",
"vplot.draw_basemap(resolution='h', ax=ax)\n",
"for c in ax.collections[-2:]:\n",
" c.set_rasterized(True)\n",
"\n",
"if False:\n",
" fig.savefig('wind-avg-power-density.pdf', transparent=True, bbox_inches='tight')"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"fig, ax = plt.subplots(figsize=(10, 6))\n",
"\n",
"ax.set_title('DE Solar average power density [GWh/a/km$^2$]')\n",
"\n",
"sh = vplot.shapes(grid_cells, solarpotential_per_grid_cell, cmap='YlOrBr', ax=ax)\n",
"vplot.shapes(network.buses.loc[onshorebuses_i[onshorede], 'region_onshore'],\n",
" outline=True, colour='black', linestyle='-', ax=ax)\n",
"vplot.shapes(network.buses.loc[offshorebuses_i[offshorede], 'region_offshore'],\n",
" outline=True, colour='black', linestyle='-', ax=ax)\n",
"fig.colorbar(sh)\n",
"\n",
"vplot.draw_basemap(resolution='h')\n",
"for c in ax.collections[-2:]:\n",
" c.set_rasterized(True)\n",
" \n",
"if False:\n",
" fig.savefig('solar-avg-power-density.pdf', transparent=True, bbox_inches='tight')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Add renewable generators to PyPSA model"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"onshorewindturbines = pd.DataFrame(dict(p_nom_extendable=True,\n",
" p_nom_max=onshore_p_nom_max,\n",
" bus=onshorebuses_i,\n",
" carrier='onwind',\n",
" marginal_cost=costs.at['onwind', 'marginal'],\n",
" capital_cost=costs.at['onwind', 'capital'],\n",
" weight=onshore_weights),\n",
" index=onshore_i)\n",
"pypsa.io.import_components_from_dataframe(network, onshorewindturbines, 'Generator')\n",
"pypsa.io.import_series_from_dataframe(\n",
" network,\n",
" onshore.sel(time=network.snapshots).transpose('time', 'name').to_pandas(),\n",
" 'Generator', 'p_max_pu'\n",
")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"offshorewindturbines = pd.DataFrame(dict(p_nom_extendable=True,\n",
" p_nom_max=offshore_p_nom_max,\n",
" bus=offshorebuses_i,\n",
" carrier='offwind',\n",
" marginal_cost=costs.at['offwind', 'marginal'],\n",
" capital_cost=costs.at['offwind', 'capital'],\n",
" weight=offshore_weights),\n",
" index=offshore_i)\n",
"pypsa.io.import_components_from_dataframe(network, offshorewindturbines, 'Generator')\n",
"pypsa.io.import_series_from_dataframe(\n",
" network,\n",
" offshore.sel(time=network.snapshots).transpose('time', 'name').to_pandas(),\n",
" 'Generator', 'p_max_pu'\n",
")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"solarpanels = pd.DataFrame(dict(p_nom_extendable=True,\n",
" p_nom_max=solar_p_nom_max,\n",
" bus=onshorebuses_i,\n",
" carrier='solar',\n",
" marginal_cost=costs.at['solar', 'marginal'],\n",
" capital_cost=costs.at['solar', 'capital'],\n",
" weight=solar_weights),\n",
" index=solar_i)\n",
"pypsa.io.import_components_from_dataframe(network, solarpanels, 'Generator')\n",
"pypsa.io.import_series_from_dataframe(\n",
" network,\n",
" solar.sel(time=network.snapshots).transpose('time', 'name').to_pandas(),\n",
" 'Generator', 'p_max_pu'\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Load"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from vresutils import load as vload"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"%%time\n",
"loads = vload.timeseries_shapes(network.buses.loc[onshorebuses_i, 'region_onshore'],\n",
" network.buses.loc[onshorebuses_i, 'country'])"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"pypsa.io.import_components_from_dataframe(network, pd.DataFrame(dict(bus=onshorebuses_i), index=onshorebuses_i), \"Load\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"pypsa.io.import_series_from_dataframe(network, loads, 'Load', 'p_set')"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Fill up the one hour\n",
"network.loads_t.p_set.ffill(inplace=True)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Export model"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"region_onshore = gpd.GeoDataFrame(dict(geometry=network.buses.region_onshore)).dropna().reset_index()\n",
"region_offshore = gpd.GeoDataFrame(dict(geometry=network.buses.region_offshore)).dropna().reset_index()\n",
"\n",
"def export_as_geojson(df, fn):\n",
" if os.path.exists(fn):\n",
" os.unlink(fn)\n",
" df.to_file(fn, driver='GeoJSON')\n",
"\n",
"export_as_geojson(region_onshore, 'results/region_onshore.json')\n",
"export_as_geojson(region_offshore, 'results/region_offshore.json')"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"network.buses['region_onshore'] = gpd.read_file('results/region_onshore.json').set_index('name')['geometry']\n",
"network.buses['region_offshore'] = gpd.read_file('results/region_offshore.json').set_index('name')['geometry']"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"network.buses.drop(['region_onshore', 'region_offshore'], axis=1, inplace=True)\n",
"network.export_to_csv_folder(\"results/02-pypsa-europe\")"
]
}
],
"metadata": {
"anaconda-cloud": {},
"git": {
"suppress_outputs": true
},
"hide_input": false,
"kernelspec": {
"display_name": "Python 2",
"language": "python",
"name": "python2"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 2
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython2",
"version": "2.7.12"
},
"nav_menu": {},
"notify_time": "5",
"toc": {
"navigate_menu": true,
"number_sections": true,
"sideBar": true,
"threshold": 6,
"toc_cell": false,
"toc_section_display": "block",
"toc_window_display": false
},
"toc_position": {
"height": "808px",
"left": "42px",
"right": "20px",
"top": "78px",
"width": "258px"
}
},
"nbformat": 4,
"nbformat_minor": 1
}