From 9aa154bd048c02fb511e78b1dc4f9524ef09251b Mon Sep 17 00:00:00 2001
From: Koen van Greevenbroek <koen.v.greevenbroek@uit.no>
Date: Thu, 25 Jan 2024 10:20:19 +0100
Subject: [PATCH] Add several map projection options

---
 config/config.default.yaml    |  1 +
 doc/configtables/plotting.csv |  1 +
 doc/release_notes.rst         |  2 ++
 scripts/plot_network.py       | 29 ++++++++++++++++++++++++-----
 4 files changed, 28 insertions(+), 5 deletions(-)

diff --git a/config/config.default.yaml b/config/config.default.yaml
index 255e3869..ef410ad6 100644
--- a/config/config.default.yaml
+++ b/config/config.default.yaml
@@ -773,6 +773,7 @@ plotting:
     color_geomap:
       ocean: white
       land: white
+  projection: "EqualEarth"
   eu_node_location:
     x: -5.5
     y: 46.
diff --git a/doc/configtables/plotting.csv b/doc/configtables/plotting.csv
index ed5d9c9f..656ed9be 100644
--- a/doc/configtables/plotting.csv
+++ b/doc/configtables/plotting.csv
@@ -1,6 +1,7 @@
 ,Unit,Values,Description
 map,,,
 -- boundaries,°,"[x1,x2,y1,y2]",Boundaries of the map plots in degrees latitude (y) and longitude (x)
+projection,--,"{EqualEarth, EuroPP, LambertAzimuthalEqualArea, LambertConformal, Orthographic}",Projection to use for maps; default is EqualEarth. LambertConformal is recommended by the European Environmental Agency.
 costs_max,bn Euro,float,Upper y-axis limit in cost bar plots.
 costs_threshold,bn Euro,float,Threshold below which technologies will not be shown in cost bar plots.
 energy_max,TWh,float,Upper y-axis limit in energy bar plots.
diff --git a/doc/release_notes.rst b/doc/release_notes.rst
index f8f958cd..0386cf49 100644
--- a/doc/release_notes.rst
+++ b/doc/release_notes.rst
@@ -35,6 +35,8 @@ Upcoming Release
 
 * Add support for the linopy ``io_api`` option; set to ``"direct"`` to increase model reading and writing performance for the highs and gurobi solvers.
 
+* Add several map projection options for plotting.
+
 
 PyPSA-Eur 0.9.0 (5th January 2024)
 ==================================
diff --git a/scripts/plot_network.py b/scripts/plot_network.py
index 33312a44..c078e412 100644
--- a/scripts/plot_network.py
+++ b/scripts/plot_network.py
@@ -170,7 +170,7 @@ def plot_map(
     line_widths = line_widths.replace(line_lower_threshold, 0)
     link_widths = link_widths.replace(line_lower_threshold, 0)
 
-    fig, ax = plt.subplots(subplot_kw={"projection": ccrs.EqualEarth()})
+    fig, ax = plt.subplots(subplot_kw={"projection": proj})
     fig.set_size_inches(7, 6)
 
     n.plot(
@@ -358,7 +358,6 @@ def plot_h2_map(network, regions):
     n.links.bus0 = n.links.bus0.str.replace(" H2", "")
     n.links.bus1 = n.links.bus1.str.replace(" H2", "")
 
-    proj = ccrs.EqualEarth()
     regions = regions.to_crs(proj.proj4_init)
 
     fig, ax = plt.subplots(figsize=(7, 6), subplot_kw={"projection": proj})
@@ -568,7 +567,7 @@ def plot_ch4_map(network):
         "biogas": "seagreen",
     }
 
-    fig, ax = plt.subplots(figsize=(7, 6), subplot_kw={"projection": ccrs.EqualEarth()})
+    fig, ax = plt.subplots(figsize=(7, 6), subplot_kw={"projection": proj})
 
     n.plot(
         bus_sizes=bus_sizes,
@@ -679,7 +678,7 @@ def plot_map_without(network):
     # Drop non-electric buses so they don't clutter the plot
     n.buses.drop(n.buses.index[n.buses.carrier != "AC"], inplace=True)
 
-    fig, ax = plt.subplots(figsize=(7, 6), subplot_kw={"projection": ccrs.EqualEarth()})
+    fig, ax = plt.subplots(figsize=(7, 6), subplot_kw={"projection": proj})
 
     # PDF has minimum width, so set these to zero
     line_lower_threshold = 200.0
@@ -993,7 +992,7 @@ def plot_map_perfect(
     link_widths[link_widths > line_upper_threshold] = line_upper_threshold
 
     for year in costs.columns:
-        fig, ax = plt.subplots(subplot_kw={"projection": ccrs.PlateCarree()})
+        fig, ax = plt.subplots(subplot_kw={"projection": proj})
         fig.set_size_inches(7, 6)
         fig.suptitle(year)
 
@@ -1082,6 +1081,26 @@ if __name__ == "__main__":
     if map_opts["boundaries"] is None:
         map_opts["boundaries"] = regions.total_bounds[[0, 2, 1, 3]] + [-1, 1, -1, 1]
 
+    proj_str = snakemake.params.plotting.get("projection", "EqualEarth")
+    central_coords = dict(central_longitude=10.0, central_latitude=50.0)
+    if proj_str == "EqualEarth":
+        # Equal area but large distortions towards the poles.
+        proj = ccrs.EqualEarth()
+    elif proj_str == "EuroPP":
+        # UTM Zone 32 projection
+        proj = ccrs.EuroPP()
+    elif proj_str == "LambertConformal":
+        # The European Environment Agency recommends using this
+        # projection for conformal pan-European mapping
+        proj = ccrs.LambertConformal(standard_parallels=(35, 65), **central_coords)
+    elif proj_str == "Orthographic":
+        proj = ccrs.Orthographic(**central_coords)
+    else:
+        logger.warning(
+            f"Plotting project {proj_str} not recognised; falling back on EqualEarth"
+        )
+        proj = ccrs.EqualEarth()
+
     if snakemake.params["foresight"] == "perfect":
         plot_map_perfect(
             n,