Japan Energy Model#
Exercise Overview#
In this exercise, we demonstrate how a workflow can be prepared using data from the Japan Energy Database.
The learning outcomes are to:
Review open-source energy data available for Japan.
Identify how raw datasets must be transformed before they can be integrated into an energy system model.
Use GeoPandas to visualise spatial data and perform basic geographic calculations.
Use pandas and Python functions to transform datasets programmatically.
Sequentially construct an energy system model using PyPSA.
Inspect and explore the resulting model structure and outputs.
References#
Ono, R., Onodera, H., Kikuchi, K. et al. (2026) ‘Japan Energy Database with hourly and municipal-scale estimates to support local energy system analysis’, Scientific Data, 13, 3. Available at: https://doi.org/10.1038/s41597-025-06294-w.
Japan Electric Power Information Center (2019) The electric power industry in Japan 2019. Tokyo: Japan Electric Power Information Center.
United Nations Office for the Coordination of Humanitarian Affairs, Regional Office for Asia and the Pacific (2019) Japan – Subnational Administrative Boundaries (administrative levels 0–2). Available at: https://data.humdata.org/dataset/cod-xa-jpn (Accessed: 6 March 2026).
lesson_number = 6
print(f'lesson{lesson_number}')
lesson6
Prepare Google Colab Environment#
#@title Connect to Google Drive {display-mode:"form"}
CONNECT_TO_DRIVE = False #@param {type:"boolean"}
import os
if CONNECT_TO_DRIVE:
from google.colab import drive
# Mount Google Drive
drive.mount('/content/drive')
# Define the desired working directory path
working_dir = '/content/drive/MyDrive/hello-pypsa'
# Create the directory if it doesn't exist
if not os.path.exists(working_dir):
os.makedirs(working_dir)
print(f"Directory '{working_dir}' created.")
else:
print(f"Directory '{working_dir}' already exists.")
# Change the current working directory
os.chdir(working_dir)
print(f"Current working directory: {os.getcwd()}")
else:
print("Not connecting to Google Drive.")
import os
#@title Install Packages {display-mode:"form"}
INSTALL_PACKAGES = False #@param {type:"boolean"}
# Check if packages have already been installed in this session to prevent re-installation
if INSTALL_PACKAGES and not os.environ.get('PYPSA_PACKAGES_INSTALLED'):
!pip install pypsa pypsa[excel] folium mapclassify cartopy
!pip install git+https://github.com/PriyeshGosai/pypsa_network_viewer.git
os.environ['PYPSA_PACKAGES_INSTALLED'] = 'true'
elif not INSTALL_PACKAGES:
print("Skipping package installation.")
else:
print("PyPSA packages are already installed for this session.")
#@title Download the file for this notebook {display-mode:"form"}
DOWNLOAD_FILE = False #@param {type:"boolean"}
if DOWNLOAD_FILE:
!gdown "https://drive.google.com/uc?id=1My8j2qRcjjhVhC5bL657oYTk7-OKDxkE"
else:
print("Skipping file download.")
Japan Model#
import pandas as pd
import geopandas as gpd
import matplotlib.pyplot as plt
import numpy as np
from unidecode import unidecode
pd.options.plotting.backend = 'plotly'
from scripts.renewables import aggregate_timeseries_by_region
from scripts.renewables import calculate_hourly_energy
print("✓ Libraries imported successfully")
japan_energy_db_file_name = r"data\JapanEnergyDatabase\JapanEnergyDatabase\JapanEnergyDatabase_2019_v1.csv"
japan_energy_db = pd.read_csv(japan_energy_db_file_name, encoding="cp932", index_col=0)
# japan_energy_db = pd.read_csv(japan_energy_db_file_name, encoding="cp932", index_col=[1, 4, 0])
japan_regions = {
"Hokkaido": {"Hokkaido"},
"Tohoku": {"Aomori", "Iwate", "Miyagi", "Akita", "Yamagata", "Fukushima"},
"Kanto": {"Tokyo", "Kanagawa", "Chiba", "Saitama", "Ibaraki", "Tochigi", "Gunma"},
"Chubu": {"Aichi", "Gifu", "Shizuoka", "Nagano", "Yamanashi", "Niigata", "Toyama", "Ishikawa", "Fukui"},
"Kansai": {"Osaka", "Kyoto", "Hyogo", "Nara", "Shiga", "Wakayama", "Mie"},
"Chugoku": {"Hiroshima", "Okayama", "Shimane", "Tottori", "Yamaguchi"},
"Shikoku": {"Ehime", "Kagawa", "Kochi", "Tokushima"},
"Kyushu": {"Fukuoka", "Saga", "Nagasaki", "Kumamoto", "Oita", "Miyazaki", "Kagoshima", "Okinawa"}
}
japan_energy_db.tail(60)
geo_data_file_name = r'data\jpn_adm_2019_shp\jpn_admbnda_adm2_2019.shp'
gdf = gpd.read_file(geo_data_file_name)
gdf
# File paths
electricity_demand_variation_file_name = r'data\JapanEnergyDatabase\Time-series\Electricity_demand_variation.csv'
offshore_wind_file_name = r'data\JapanEnergyDatabase\Time-series\Offshore_wind.csv'
onshore_wind_file_name = r'data\JapanEnergyDatabase\Time-series\Onshore_wind.csv'
photovoltaic_file_name = r'data\JapanEnergyDatabase\Time-series\Photovoltaic.csv'
runofriver_file_name = r'data\JapanEnergyDatabase\Time-series\Runofriver.csv'
energy_capacity = pd.read_csv('data/energy_capacity.csv', index_col=[0, 1]) # energy_capacity.xs('Total', level=1)['2016']
# Load all data files
# japan_energy_db = pd.read_csv(japan_energy_db_file_name, encoding="cp932", index_col=0)
# Load timeseries files with proper datetime parsing
electricity_demand_variation = pd.read_csv(electricity_demand_variation_file_name, index_col=0)
electricity_demand_variation.index = pd.to_datetime(electricity_demand_variation.index, format='mixed', dayfirst=True)
electricity_demand_variation.index.name = None
offshore_wind = pd.read_csv(offshore_wind_file_name, index_col=0)
offshore_wind.index = pd.to_datetime(offshore_wind.index, format='mixed', dayfirst=True)
offshore_wind.index.name = None
onshore_wind = pd.read_csv(onshore_wind_file_name, index_col=0)
onshore_wind.index = pd.to_datetime(onshore_wind.index, format='mixed', dayfirst=True)
onshore_wind.index.name = None
photovoltaic = pd.read_csv(photovoltaic_file_name, index_col=0)
photovoltaic.index = pd.to_datetime(photovoltaic.index, format='mixed', dayfirst=True)
photovoltaic.index.name = None
runofriver = pd.read_csv(runofriver_file_name, index_col=0)
runofriver.index = pd.to_datetime(runofriver.index, format='mixed', dayfirst=True)
runofriver.index.name = None
print(f"✓ Data files loaded successfully")
print(f" - japan_energy_db: {japan_energy_db.shape}")
print(f" - Timeseries data shape: {onshore_wind.shape}")
print(f" - Timeseries index type: {type(onshore_wind.index)}")
print(f" - First timestamp: {onshore_wind.index[0]}")
print(f" - Last timestamp: {onshore_wind.index[-1]}")
electricity_demand_variation.head()
japan_energy_db
gdf
2. Translation and Mapping Setup#
# Normalization function for consistent string matching
def normalize_name(s):
"""Normalize names by removing spaces, accents, and converting to lowercase"""
if isinstance(s, str):
s = s.replace("\xa0", " ") # remove non-breaking spaces
s = unidecode(s) # remove accents/macrons
s = s.strip() # remove extra spaces
s = s.lower() # standardize case
return s
# Apply normalization
gdf["prefecture_norm"] = gdf["ADM1_EN"].apply(normalize_name)
print("✓ Normalization function defined and applied")
# Build japan_translation dictionary from GDF data (NOT hardcoded!)
# Only 8 region translations need to be hardcoded
country_translations = gdf[['ADM0_JA', 'ADM0_EN']].drop_duplicates()
prefecture_translations = gdf[['ADM1_JA', 'ADM1_EN']].drop_duplicates()
japan_translation = {
"country": dict(zip(country_translations['ADM0_JA'], country_translations['ADM0_EN'])),
"regions": {
# Only these 8 region names are hardcoded (not in gdf)
"東北": "Tohoku",
"関東": "Kanto",
"北陸": "Hokuriku",
"中部": "Chubu",
"関西": "Kansai",
"中国": "Chugoku",
"四国": "Shikoku",
"九州": "Kyushu"
},
"prefectures": dict(zip(prefecture_translations['ADM1_JA'], prefecture_translations['ADM1_EN']))
}
# Flatten for easy mapping
japan_translation_flat = {}
for category in japan_translation.values():
if isinstance(category, dict):
japan_translation_flat.update(category)
print(f"✓ Translation dictionary built from GDF data:")
print(f" - Country: {len(japan_translation['country'])} (from gdf)")
print(f" - Regions: {len(japan_translation['regions'])} (hardcoded)")
print(f" - Prefectures: {len(japan_translation['prefectures'])} (from gdf)")
print(f" - Total hardcoded items: {len(japan_translation['regions'])} (was 103!)")
# PCODE-based region mapping (replaces 47 hardcoded prefecture→region mappings)
def get_region_from_pcode(pcode):
"""
Map prefecture PCODE to region using Japanese administrative numbering.
This eliminates the need for 47 hardcoded prefecture→region mappings!
"""
pcode_num = int(pcode.replace('JP', ''))
if pcode_num == 1:
return 'Hokkaido'
elif 2 <= pcode_num <= 7:
return 'Tohoku'
elif 8 <= pcode_num <= 14:
return 'Kanto'
elif 15 <= pcode_num <= 18:
return 'Hokuriku'
elif 19 <= pcode_num <= 23:
return 'Chubu'
elif 24 <= pcode_num <= 30:
return 'Kansai'
elif 31 <= pcode_num <= 35:
return 'Chugoku'
elif 36 <= pcode_num <= 39:
return 'Shikoku'
elif 40 <= pcode_num <= 47:
return 'Kyushu'
else:
return 'Unknown'
# Apply region mapping
gdf['region'] = gdf['ADM1_PCODE'].apply(get_region_from_pcode)
print(f"✓ PCODE-based region mapping applied:")
print(f" - Regions: {sorted(gdf['region'].unique())}")
print(f" - All {len(gdf)} entries mapped successfully")
# Apply translations to japan_energy_db
japan_energy_db['pref_en'] = japan_energy_db['pref'].map(japan_translation_flat)
# Categorize entries (country, region, or prefecture)
def get_category(japanese_name):
if pd.isna(japanese_name):
return None
if japanese_name in japan_translation['country']:
return 'country'
elif japanese_name in japan_translation['regions']:
return 'region'
elif japanese_name in japan_translation['prefectures']:
return 'prefecture'
else:
return 'unknown'
japan_energy_db['category'] = japan_energy_db['pref'].apply(get_category)
3. Region-Level Analysis#
# Create region-level geometries by dissolving prefecture boundaries
regions_gdf = gdf.dissolve(by='region', as_index=False)
regions_gdf['centroid'] = regions_gdf.geometry.centroid
# Create region adjacency and distance matrix
def are_regions_adjacent(geom1, geom2):
"""Check if two geometries share a border"""
return geom1.touches(geom2) or geom1.intersects(geom2)
def calculate_centroid_distance(centroid1, centroid2):
"""Calculate distance between centroids in kilometers"""
distance = centroid1.distance(centroid2)
return distance * 111 # Approximate conversion from degrees to km
# Special connections (e.g., Seikan Tunnel)
special_connections = [('Hokkaido', 'Tohoku'), ('Shikoku', 'Kansai'), ('Kyushu','Chugoku')]
# Build distance matrix
region_names = sorted(regions_gdf['region'].tolist())
distance_matrix_data = {}
for region1 in region_names:
row_data = {}
geom1 = regions_gdf[regions_gdf['region'] == region1].iloc[0].geometry
centroid1 = regions_gdf[regions_gdf['region'] == region1].iloc[0].centroid
for region2 in region_names:
geom2 = regions_gdf[regions_gdf['region'] == region2].iloc[0].geometry
centroid2 = regions_gdf[regions_gdf['region'] == region2].iloc[0].centroid
if region1 == region2:
row_data[region2] = 0
elif are_regions_adjacent(geom1, geom2) or \
(region1, region2) in special_connections or \
(region2, region1) in special_connections:
distance_km = calculate_centroid_distance(centroid1, centroid2)
row_data[region2] = round(distance_km, 2)
else:
row_data[region2] = 'X'
distance_matrix_data[region1] = row_data
distance_matrix_df = pd.DataFrame(distance_matrix_data).T
4. Visualizations#
# Visualize regions on map
fig, ax = plt.subplots(1, 1, figsize=(14, 12))
regions_gdf.plot(ax=ax, edgecolor='black', linewidth=1.5, alpha=0.6,
column='region', categorical=True, legend=True, cmap='Set3')
# Add region labels
for idx, row in regions_gdf.iterrows():
centroid = row.centroid
ax.annotate(text=row['region'], xy=(centroid.x, centroid.y),
fontsize=11, ha='center', va='center', fontweight='bold',
bbox=dict(boxstyle='round,pad=0.5', facecolor='white', alpha=0.8))
ax.set_title('Japan Regions Map', fontsize=16, fontweight='bold', pad=20)
ax.axis('off')
plt.tight_layout()
plt.show()
# Visualize distance matrix as heatmap
fig, ax = plt.subplots(figsize=(12, 10))
# Convert to numeric for visualization
distance_matrix_numeric = distance_matrix_df.copy()
for i in range(len(distance_matrix_numeric)):
for j in range(len(distance_matrix_numeric.columns)):
val = distance_matrix_numeric.iloc[i, j]
if val == 'X':
distance_matrix_numeric.iloc[i, j] = np.nan
else:
distance_matrix_numeric.iloc[i, j] = float(val)
distance_matrix_numeric = distance_matrix_numeric.astype(float)
# Create heatmap
im = ax.imshow(distance_matrix_numeric, cmap='YlOrRd', aspect='auto')
ax.set_xticks(np.arange(len(region_names)))
ax.set_yticks(np.arange(len(region_names)))
ax.set_xticklabels(region_names, rotation=45, ha='right')
ax.set_yticklabels(region_names)
# Add text annotations
for i in range(len(region_names)):
for j in range(len(region_names)):
val = distance_matrix_df.iloc[i, j]
if val == 0:
text = '-'
color = 'black'
elif val == 'X':
text = 'X'
color = 'gray'
else:
text = f'{val}'
color = 'white' if val > 200 else 'black'
ax.text(j, i, text, ha='center', va='center', color=color, fontsize=8)
ax.set_title('Region Adjacency & Centroid Distance Matrix (km)\n"X" = Not Adjacent',
fontsize=14, fontweight='bold', pad=20)
plt.colorbar(im, ax=ax, label='Distance (km)')
plt.tight_layout()
plt.show()
5. Data Processing#
# Load all data files
common_headings = ['pref',
'city',
'kana',
'roma']
electric_demand = ['ele_ind',
'ele_com',
'ele_res',
'ele_tra']
potential_re = ['potential_onwind',
'potential_offwind',
'potential_solar_roof',
'potential_solar_util',
'potential_solar',
'potential_river']
potential_conv = ['potential_geotherm',
'potential_biomass',
'potential_waste']
exist_re = ['exist_onwind',
'exist_offwind',
'exist_solar',
'exist_solar_roof',
'exist_solar_util',
'exist_river']
exist_conv = ['exist_geotherm',
'exist_biomass',
'exist_waste']
# Dictionary mapping technology columns to their timeseries DataFrames
# Note: potential_solar and potential_river don't have matching timeseries in current data
potential_re_df = japan_energy_db.loc[:, potential_re].copy()
# Convert index to string to match offshore_wind columns
potential_re_df.index = potential_re_df.index.astype(str)
re_technology_mapping = {
'potential_onwind': onshore_wind,
'potential_offwind': offshore_wind,
'potential_solar': photovoltaic, # Using photovoltaic for all solar
'potential_solar_roof': photovoltaic,
'potential_solar_util': photovoltaic,
'potential_river': runofriver
}
# Calculate hourly energy for all renewable technologies
print("=== Calculating hourly energy for all RE technologies ===\n")
re_hourly_energy = {}
for tech_col, timeseries_df in re_technology_mapping.items():
print(f"\n--- {tech_col} ---")
try:
energy_df = calculate_hourly_energy(
capacity_df=potential_re_df,
timeseries_df=timeseries_df,
technology_column=tech_col
)
re_hourly_energy[tech_col] = energy_df
except Exception as e:
print(f"✗ Error processing {tech_col}: {e}")
print(f"\n=== Summary ===")
print(f"Successfully processed {len(re_hourly_energy)} technologies:")
for tech_name, energy_df in re_hourly_energy.items():
print(f" {tech_name}: {energy_df.shape[1]} nodes, {energy_df.shape[0]} timesteps")
# Prepare existing RE capacity DataFrame
exist_re_df = japan_energy_db.loc[:, exist_re].copy()
exist_re_df.index = exist_re_df.index.astype(str)
# Dictionary mapping existing technology columns to their timeseries DataFrames
exist_re_technology_mapping = {
'exist_onwind': onshore_wind,
'exist_offwind': offshore_wind,
'exist_solar': photovoltaic,
'exist_solar_roof': photovoltaic,
'exist_solar_util': photovoltaic,
'exist_river': runofriver
}
# Calculate hourly energy for all existing renewable technologies
print("=== Calculating hourly energy for existing RE technologies ===\n")
exist_re_hourly_energy = {}
for tech_col, timeseries_df in exist_re_technology_mapping.items():
print(f"\n--- {tech_col} ---")
try:
energy_df = calculate_hourly_energy(
capacity_df=exist_re_df,
timeseries_df=timeseries_df,
technology_column=tech_col
)
exist_re_hourly_energy[tech_col] = energy_df
except Exception as e:
print(f"✗ Error processing {tech_col}: {e}")
print(f"\n=== Summary ===")
print(f"Successfully processed {len(exist_re_hourly_energy)} existing RE technologies:")
for tech_name, energy_df in exist_re_hourly_energy.items():
print(f" {tech_name}: {energy_df.shape[1]} nodes, {energy_df.shape[0]} timesteps")
from scripts.renewables import prepare_capacity_dataframes
# Map technologies to their timeseries
timeseries_mapping = {
'potential_onwind': onshore_wind,
'potential_offwind': offshore_wind,
'potential_solar': photovoltaic,
'potential_solar_roof': photovoltaic,
'potential_solar_util': photovoltaic,
'potential_river': runofriver,
'exist_onwind': onshore_wind,
'exist_offwind': offshore_wind,
'exist_solar': photovoltaic,
'exist_solar_roof': photovoltaic,
'exist_solar_util': photovoltaic,
'exist_river': runofriver
}
# Combine potential and existing into one dataframe
all_re_df = pd.concat([potential_re_df, exist_re_df], axis=1)
# Convert all TJ/year to MW
p_nom_all = prepare_capacity_dataframes(all_re_df, timeseries_mapping)
# Split into p_nom_max and p_nom
p_nom_max = p_nom_all[potential_re]
p_nom = p_nom_all[exist_re]
print("p_nom_max shape:", p_nom_max.shape)
print("p_nom shape:", p_nom.shape)
print("\nSample p_nom_max:")
print(p_nom_max.head())
print("\nSample p_nom:")
print(p_nom.head())
# Verify: Check node 11002
code = '11002'
print(f"Node {code} verification:")
print(f" p_nom_max (onwind): {p_nom_max.loc[code, 'potential_onwind']:.2f} MW")
print(f" p_nom (onwind): {p_nom.loc[code, 'exist_onwind']:.2f} MW")
print(f" Different? {p_nom_max.loc[code, 'potential_onwind'] != p_nom.loc[code, 'exist_onwind']}")
Aggregate Capacities by Region#
from scripts.renewables import create_code_to_region_mapping, aggregate_by_region
# Create code to region mapping
code_to_region = create_code_to_region_mapping(gdf)
# Aggregate to regional level
p_nom_max_regional = aggregate_by_region(p_nom_max, code_to_region)
p_nom_regional = aggregate_by_region(p_nom, code_to_region)
print("Regional p_nom_max:")
print(p_nom_max_regional)
print("\nRegional p_nom:")
print(p_nom_regional)
# Comparison: Regional potential vs existing
regional_comparison = pd.DataFrame({
'Region': p_nom_max_regional.index,
'Potential (MW)': p_nom_max_regional.sum(axis=1),
'Existing (MW)': p_nom_regional.sum(axis=1)
})
regional_comparison['Utilization (%)'] = (
regional_comparison['Existing (MW)'] / regional_comparison['Potential (MW)'] * 100
).round(2)
print(regional_comparison.to_string(index=False))
Regional Average Renewable Profiles#
from scripts.renewables import aggregate_timeseries_by_region
# Calculate regional average profiles for each renewable technology
onshore_wind_regional = aggregate_timeseries_by_region(onshore_wind, code_to_region)
offshore_wind_regional = aggregate_timeseries_by_region(offshore_wind, code_to_region)
photovoltaic_regional = aggregate_timeseries_by_region(photovoltaic, code_to_region)
runofriver_regional = aggregate_timeseries_by_region(runofriver, code_to_region)
print("Regional profiles created:")
print(f" onshore_wind_regional: {onshore_wind_regional.shape}")
print(f" offshore_wind_regional: {offshore_wind_regional.shape}")
print(f" photovoltaic_regional: {photovoltaic_regional.shape}")
print(f" runofriver_regional: {runofriver_regional.shape}")
print("\nSample onshore_wind_regional (first 5 hours):")
print(onshore_wind_regional.head())
Demand data#
sum_series = japan_energy_db.loc[:, electric_demand].sum(axis=1)
sum_series.index = sum_series.index.astype(str)
demand_df = electricity_demand_variation.multiply(
sum_series[electricity_demand_variation.columns], axis='columns'
)
# Aggregate demand by region using the same approach as renewable timeseries
demand_regional = aggregate_timeseries_by_region(demand_df, code_to_region)
print("Regional demand created:")
print(f" demand_regional: {demand_regional.shape}")
print(f" Index type: {type(demand_regional.index)}")
print(f"\nSample demand_regional (first 5 hours):")
print(demand_regional.head())
print(f"\nTotal annual demand by region (GWh):")
print((demand_regional.sum() / 1000).round(2))
6. PyPSA Network Setup#
import pypsa
# Initialize network
network = pypsa.Network()
# Get list of regions
regions = p_nom_max_regional.index.tolist()
print(f"Creating network with {len(regions)} regions: {regions}")
network.add("Carrier", "onwind", color="#4EA8DE") # steel blue
network.add("Carrier", "offwind", color="#1B6CA8") # dark blue
network.add("Carrier", "solar_roof", color="#FFD700") # gold
network.add("Carrier", "solar_util", color="#FFA500") # orange
network.add("Carrier", "solar", color="#FF8C00") # dark orange
network.add("Carrier", "river", color="#2ECC71") # green
network.add("Carrier", "geotherm", color="#E74C3C") # red
network.add("Carrier", "biomass", color="#8B6914") # brown
network.add("Carrier", "waste", color="#7D6608") # olive
# Set network snapshots to simple hour indices (0-8759)
network.set_snapshots(range(8760))
print(f"Network snapshots set to {len(network.snapshots)} hours")
print(f"Snapshots: hour {network.snapshots[0]} to hour {network.snapshots[-1]}")
network.snapshots
# Loop 1: Create buses for each region
for region in regions:
# Get centroid coordinates for this region
region_centroid = regions_gdf[regions_gdf['region'] == region].iloc[0].centroid
network.add("Bus", region,
x=region_centroid.x,
y=region_centroid.y)
# Add bus options here
print(f"Added {len(network.buses)} buses")
print("\nBus locations:")
print(network.buses[['x', 'y']])
distance_matrix_df
transmission_cost_factor = 1
# Loop 2: Create bi-directional links between connected regions
# Get connected regions from distance matrix (where value is not 'X')
for region1 in regions:
for region2 in regions:
if region1 < region2: # Avoid duplicates
distance = distance_matrix_df.loc[region1, region2]
if distance != 'X' and distance > 0:
# Forward link
network.add(
"Link",
f"{region1}-{region2}",
bus0=region1,
bus1=region2,
p_min_pu=-1,
p_nom_max = 5,
length = distance,
capital_cost = distance * transmission_cost_factor
# Add link options here
)
print(f"Added {len(network.links)} links")
# Prepare renewable timeseries mapping
renewable_timeseries = {
'onwind': onshore_wind_regional,
'offwind': offshore_wind_regional,
'solar_roof': photovoltaic_regional,
'solar_util': photovoltaic_regional,
'solar': photovoltaic_regional,
}
# Prepare capacity mapping (potential)
renewable_capacity_max = {
'onwind': 'potential_onwind',
'offwind': 'potential_offwind',
'solar_roof': 'potential_solar_roof',
'solar_util': 'potential_solar_util',
'solar': 'potential_solar',
}
# Prepare capacity mapping (existing)
renewable_capacity_nom = {
'onwind': 'exist_onwind',
'offwind': 'exist_offwind',
'solar_roof': 'exist_solar_roof',
'solar_util': 'exist_solar_util',
'solar': 'exist_solar',
}
print("Renewable timeseries prepared")
# Loop 3a: Add renewable generators (onwind, offwind, solar_roof, solar_util, solar)
for region in regions:
for tech_name, timeseries_df in renewable_timeseries.items():
capacity_col = renewable_capacity_max[tech_name]
p_nom_max_value = p_nom_max_regional.loc[region, capacity_col]
if pd.notna(p_nom_max_value) and p_nom_max_value > 0:
network.add(
"Generator",
f"{region}-{tech_name}",
bus=region,
p_max_pu=timeseries_df[region],
p_nom_max=p_nom_max_value,
# Add generator options here
)
print(f"Added {len(network.generators)} renewable generators")
# Loop 3b: Add conventional generators (geotherm, biomass, waste)
conventional_techs = ['geotherm', 'biomass', 'waste']
for region in regions:
for tech_name in conventional_techs:
capacity_col = f'potential_{tech_name}'
# Check if column exists in p_nom_max_regional
if capacity_col in p_nom_max_regional.columns:
p_nom_max_value = p_nom_max_regional.loc[region, capacity_col]
if pd.notna(p_nom_max_value) and p_nom_max_value > 0:
network.add(
"Generator",
f"{region}-{tech_name}",
bus=region,
p_nom_max=p_nom_max_value,
# Add generator options here
)
print(f"Total generators: {len(network.generators)}")
# Loop 4: Add storage units (river)
for region in regions:
capacity_col = 'potential_river'
p_nom_max_value = p_nom_max_regional.loc[region, capacity_col]
if pd.notna(p_nom_max_value) and p_nom_max_value > 0:
network.add(
"StorageUnit",
f"{region}-river",
bus=region,
p_nom_max=p_nom_max_value,
p_max_pu=runofriver_regional[region],
# Add storage unit options here
)
print(f"Added {len(network.storage_units)} storage units")
# Loop 5: Add loads
for region in regions:
network.add(
"Load",
f"{region}-load",
bus=region,
p_set=demand_regional[region]
# Add load options here
)
print(f"Added {len(network.loads)} loads")
print("\nLoad statistics:")
print(network.loads_t.p_set.sum().round(2))
# Summary
print("\n=== Network Summary ===")
print(f"Buses: {len(network.buses)}")
print(f"Links: {len(network.links)}")
print(f"Generators: {len(network.generators)}")
print(f"Storage Units: {len(network.storage_units)}")
print(f"Loads: {len(network.loads)}")
network.links