fix subsidy module

CHANGELOG.md

@@ -27,6 +27,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - Global Energy Monitor data integration
 - Technology configuration system
 - Network plotting with customizable styles
+- Policies (subsidies) and differentiated fuel costs
 
 ### Changed
 - Improved documentation structure with tutorials and reference guides
@@ -67,7 +68,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ## Version History Notes
 
-PyPSA-China (PIK) is adapted from the Zhou et al. version, which was originally developed by Hailiang Liu et al. This changelog tracks changes from version 1.0.0 onwards in the PIK implementation.
+PyPSA-China (PIK) is based on the paper by Zhou et al, which extends a version original developed by Hailiang Liu et al. This changelog tracks changes from version 1.0.0 onwards in the PIK implementation.
 
 For detailed information about specific changes, see the [commit history](https://github.com/pik-piam/PyPSA-China-PIK/commits/main) on GitHub.
 

docs/configuration.md

@@ -542,18 +542,33 @@ nodes:
 - **`splits`**: Custom groupings of admin level 2 regions within provinces
 
 ## Fuel Subsidies
-
-```yaml
-subsidies:
-  enabled: false
-  gas:
-    Guangdong: -10
-    Jiangsu: -10
-    Zhejiang: -10
-    Beijing: -11
-    Tianjin: -11
-    Shanghai: -11
-```
+Fuel subsidies can be speficied for all years or per year
+
+    === "All years"
+      ```yaml
+      subsidies:
+        enabled: false
+        gas:
+          Guangdong: -10
+          Jiangsu: -10
+          Zhejiang: -10
+          Beijing: -11
+          Tianjin: -11
+          Shanghai: -11
+      ```
+    === "Year by year"
+      ```yaml
+      subsidies:
+        enabled: false
+        gas:
+          2020:
+            Guangdong: -10
+            Jiangsu: -10
+            Zhejiang: -10
+            Beijing: -11
+            Tianjin: -11
+            Shanghai: -11
+      ```
 
 Provincial fuel subsidies configuration:
 - **`enabled`**: Enable/disable fuel subsidy system

examples/historical.yml

@@ -0,0 +1,34 @@
+# A Configuration to reproduce historical mix for 2020 and 2025
+# run with `snakemake --configfile=examples/historical.yml`
+
+run:
+  name: "reproduce_historical_load"
+foresight: "overnight"
+
+scenario: 
+  co2_pathway: ["exp175default"] # co2_scenarios that will be used
+  topology: "current+FCG" # "current" or "FCG" or "current+FCG" or "current+Neighbor"
+  planning_horizons:
+    - 2020
+    - 2025
+
+subsidies:
+  enabled: True  # Set to false to disable fuel subsidies
+  # Year-dependent format: subsidies.fuel_type -> year -> province -> value (EUR/MWh)
+  # Only negative values allowed (subsidies reduce marginal cost)
+  # Gas favoured over coal in urban areas due to PM2.5 concerns
+  gas:
+    2020:
+      Guangdong: -10.16
+      Jiangsu:  -10.16
+      Zhejiang: -10.15
+      Beijing: -12.5  
+      Tianjin: -12.5
+      Shanghai: -12.5
+      Xinjiang: -10
+  # location dependent fuel prices (cheaper in Shaanxi and Inner Mongolia for example)
+  coal:
+    2020:
+      Xinjiang: -4
+      InnerMongolia: -4.5
+      Hebei: -4

workflow/__init__.py

@@ -1,2 +1,3 @@
+"""Track version"""
 # pypsa-China PIK editions
 __version__ = "1.3.2"

workflow/scripts/__init__.py

@@ -1 +1,2 @@
+"""init file for scope"""
 # for make the docs

workflow/scripts/_plot_utilities.py

@@ -14,25 +14,25 @@
 
 def validate_hex_colors(tech_colors: dict[str, str], fill_color = "#999999") -> dict[str, str]:
     """Validate and standardize hex color codes in technology color mappings.
-    
+
     Ensures all color codes in the technology colors dictionary are valid hexadecimal
     color codes. Invalid or malformed colors are replaced with a default gray color.
-    
+
     Args:
         tech_colors (Dict[str, str]): Dictionary mapping technology names to color codes. Expected
             format is {'tech_name': '#RRGGBB'} or {'tech_name': '#RGB'}.
         fill_color (str, optional): Default color to use for invalid entries. Defaults to '#999999'.
-            
+
     Returns:
         dict[str,str] with validated hex color codes. All valid colors are converted
         to lowercase, while invalid colors are replaced with '#999999' (gray).
-        
+
     Example:
         >>> colors = {'solar': '#FFD700', 'wind': 'invalid', 'coal': '#8B4513'}
         >>> validated = validate_hex_colors(colors)
         >>> print(validated)
         {'solar': '#ffd700', 'wind': '#999999', 'coal': '#8b4513'}
-        
+
     Note:
         Accepts both 3-digit (#RGB) and 6-digit (#RRGGBB) hex color formats.
         All valid colors are standardized to lowercase.
@@ -463,16 +463,16 @@ def annotate_heatmap(
 
 def setup_plot_export_hook(plot_accessor_class, export_dir="plot_exports", verbose=True):
     """Setup a monkey patch to auto-export data to CSV whenever pandas plots are created.
-    
+
     Args:
         plot_accessor_class: The PlotAccessor class to patch (e.g., pandas.plotting.PlotAccessor).
         export_dir (str, optional): Directory where CSV exports will be saved. 
             Defaults to "plot_exports".
         verbose (bool, optional): Whether to print export messages. Defaults to True.
-    
+
     Returns:
         callable: Function to remove the patch and restore original behavior.
-    
+
     Example:
         >>> from pandas.plotting import PlotAccessor
         >>> remove_hook = setup_plot_export_hook(PlotAccessor)
@@ -481,15 +481,16 @@ def setup_plot_export_hook(plot_accessor_class, export_dir="plot_exports", verbo
     """
     import os
     import time
+
     import pandas as pd
-    
+
     # Create export directory
     os.makedirs(export_dir, exist_ok=True)
-    
+
     # Store original __call__ if not already stored
     if not hasattr(plot_accessor_class, '_original_call'):
         plot_accessor_class._original_call = plot_accessor_class.__call__
-    
+
     def patched_plot_call(self, *args, **kwargs):
         """Patched __call__ method for PlotAccessor to export data before plotting."""
         # Create timestamped filename
@@ -498,22 +499,22 @@ def patched_plot_call(self, *args, **kwargs):
         else:
             ts = time.strftime("%Y%m%d_%H%M%S")
             fname = os.path.join(export_dir, f"plot_export_{ts}.csv")
-        
+
         # Export the data
         if isinstance(self._parent, pd.Series):
             self._parent.to_frame().to_csv(fname, index=True)
         else:
             self._parent.to_csv(fname, index=True)
-        
+
         if verbose:
             print(f"[pandas-plot-hook] Exported plotted data to {fname}")
-        
+
         # Call the original __call__ method
         return self._original_call(*args, **kwargs)
-    
+
     # Apply the patch
     plot_accessor_class.__call__ = patched_plot_call
-    
+
     # Return function to remove the patch
     def remove_hook():
         """Remove the plot export hook and restore original behavior."""
@@ -522,7 +523,7 @@ def remove_hook():
             delattr(plot_accessor_class, '_original_call')
             if verbose:
                 print("[pandas-plot-hook] Hook removed, original behavior restored.")
-    
+
     return remove_hook
 
 

workflow/scripts/_pypsa_helpers.py

@@ -8,8 +8,8 @@
 import pandas as pd
 import pypsa
 import pytz
-
 from constants import PROV_NAMES
+
 # get root logger
 logger = logging.getLogger()
 

workflow/scripts/add_sectors.py

@@ -11,7 +11,7 @@
 
 def add_carrier_if_missing(n: pypsa.Network, carrier_name: str):
     """Add a carrier to the network if it doesn't already exist.
-    
+
     Args:
         n (pypsa.Network): PyPSA network to modify.
         carrier_name (str): Name of the carrier to add.
@@ -113,7 +113,7 @@ def attach_simple_ev(
 
     transport_cfg = ev_cfg.get("transport", {})
     logger.info("Transport configuration: %s", transport_cfg)
-    
+
     passenger_cfg = transport_cfg.get("passenger_bev", {})
     if passenger_cfg.get("enable", False):
         charging = pd.read_csv(

workflow/scripts/build_population.py

@@ -14,7 +14,7 @@
 
 def load_pop_csv(csv_path: os.PathLike) -> pd.DataFrame:
     """Load the national bureau of statistics of China population.
-    
+
     Supports both formats:
     - Yearbook format (2.5 pop at year end by Region)
     - Historical data format with comment lines
@@ -24,7 +24,7 @@ def load_pop_csv(csv_path: os.PathLike) -> pd.DataFrame:
 
     Returns:
         pd.DataFrame: The population for constants.POP_YEAR by province
-        
+
     Raises:
         ValueError: If the province names do not match expected names
     """

workflow/scripts/determine_availability_matrix.py

@@ -27,7 +27,6 @@
 import numpy as np
 from _helpers import configure_logging, mock_snakemake
 from constants import OFFSHORE_WIND_NODES, PROV_NAMES
-from pandas import concat
 from readers_geospatial import read_offshore_province_shapes, read_province_shapes
 
 logger = logging.getLogger(__name__)

workflow/scripts/functions.py

@@ -19,28 +19,28 @@
 # polynomial centroid for plotting
 def get_poly_center(poly: Polygon):
     """Get the geographic centroid of a polygon geometry.
-    
+
     Extracts the centroid coordinates from a polygon object, typically used
     for plotting and spatial analysis in geographic applications.
-    
+
     Args:
         poly (Polygon): A (shapely) polygon geometry object with a 
             centroid attribute that has x and y coordinate arrays.
-            
+
     Returns:
         tuple: A tuple containing (x, y) coordinates of the polygon centroid
             as floating point numbers.
-            
+
     Example:
         >>> from shapely.geometry import Polygon
         >>> polygon = Polygon([(0, 0), (1, 0), (1, 1), (0, 1)])
         >>> center = get_poly_center(polygon)
         >>> print(center)
         (0.5, 0.5)
-        
+
     Note:
-        This function assumes the polygon object has a centroid attribute
-        with xy arrays containing coordinate data.
+    This function assumes the polygon object has a centroid attribute
+    with xy arrays containing coordinate data.
     """
     return (poly.centroid.xy[0][0], poly.centroid.xy[1][0])
 
@@ -69,26 +69,26 @@ def cartesian(s1: pd.Series, s2: pd.Series) -> pd.DataFrame:
 
 def haversine(p1, p2) -> float:
     """Calculate the great circle distance between two points on Earth.
-    
+
     Uses the Haversine formula to compute the shortest distance over the Earth's
     surface between two points specified in decimal degrees latitude and longitude.
     This is useful for calculating distances between geographic locations.
-    
+
     Args:
         p1 (shapely.Point): location 1 in decimal deg
         p2 (shapely.Point): location 2 in decimal deg
 
     Returns:
         float: Great circle distance between the two points in kilometers.
-        
+
     Example:
         >>> from shapely.geometry import Point
         >>> beijing = Point(116.4074, 39.9042)  # longitude, latitude
         >>> shanghai = Point(121.4737, 31.2304)
         >>> distance = haversine(beijing, shanghai)
         >>> print(f"Distance: {distance:.1f} km")
         Distance: 1067.1 km
-        
+
     Note:
         The function assumes the Earth is a perfect sphere with radius 6371 km.
     """
@@ -112,7 +112,7 @@ def area_from_lon_lat_poly(geometry: Polygon):
 
     Args:
         geometry (Polygon): Polygon geometry in lon-lat coordinates.
-    
+
     Returns:
         float: Area of the polygon in m^2.
     """

workflow/scripts/plot_network.py

@@ -732,4 +732,4 @@ def calc_plot_width(row, carrier="AC"):
         save_path=snakemake.output.cost_map,
     )
 
-    logger.info("Network successfully plotted")
+    logger.info("Network successfully plotted")

workflow/scripts/plot_time_series.py

@@ -267,10 +267,9 @@ def plot_residual_load_duration_curve(
         )
         .groupby(level=1)
         .sum()
-        .loc[vre_techs]
-        .sum()
     )
-
+    tech_filter = [t for t in vre_techs if t in vre_supply.index]
+    vre_supply = vre_supply.loc[tech_filter].sum()
     residual = (load - vre_supply).sort_values(ascending=False) / PLOT_CAP_UNITS
     residual.reset_index(drop=True).plot(ax=ax, lw=3)
     ax.set_ylabel(f"Residual Load [{PLOT_CAP_LABEL}]")
@@ -529,7 +528,7 @@ def plot_vre_timemap(
             # co2_pathway="SSP2-PkBudg1000-CHA-pypsaelh2",
             heating_demand="positive",
             # configfiles=["resources/tmp/remind_coupled_cg.yaml"],
-            planning_horizons="2050",
+            planning_horizons="2025",
             winter_day1="12-10 21:00",  # mm-dd HH:MM
             winter_day2="12-17 12:00",  # mm-dd HH:MM
             spring_day1="03-31 21:00",  # mm-dd HH:MM