From 88dd5ef9d86886d04d57c1b83a8036161c596ea4 Mon Sep 17 00:00:00 2001 From: Amos Schledorn Date: Mon, 22 Dec 2025 16:33:29 +0100 Subject: [PATCH 01/22] chor: integrate lake heat into existing structure --- data/versions.csv | 2 + rules/build_sector.smk | 35 ++ rules/retrieve.smk | 30 ++ .../lake_water_heat_approximator.py | 114 +++++ .../river_water_heat_approximator.py | 3 +- .../build_lake_water_heat_potential.py | 444 ++++++++++++++++++ 6 files changed, 627 insertions(+), 1 deletion(-) create mode 100644 scripts/build_surface_water_heat_potentials/approximators/lake_water_heat_approximator.py create mode 100644 scripts/build_surface_water_heat_potentials/build_lake_water_heat_potential.py diff --git a/data/versions.csv b/data/versions.csv index e18f2ea663..a72c678ade 100644 --- a/data/versions.csv +++ b/data/versions.csv @@ -99,3 +99,5 @@ "geothermal_heat_utilisation_potentials","archive","341.5","['latest', 'supported']","https://zenodo.org/records/17207640/files/Results_DH_Matching_Cluster.xlsx","" "jrc_ardeco","primary","2021","['latest', 'supported']","https://territorial.ec.europa.eu/ardeco-api-v2/rest/export/","" "jrc_ardeco","archive","2021","['latest', 'supported']","https://zenodo.org/records/17249457/files","" +"lake_data","primary","v10","['latest', 'supported']",""https://data.hydrosheds.org/file/hydrolakes/HydroLAKES_polys_v10.gdb.zip"","" +"lake_data","archive","v10","['latest', 'supported']",""https://data.hydrosheds.org/file/hydrolakes/HydroLAKES_polys_v10.gdb.zip"","" \ No newline at end of file diff --git a/rules/build_sector.smk b/rules/build_sector.smk index bd239b2584..c519c121df 100755 --- a/rules/build_sector.smk +++ b/rules/build_sector.smk @@ -496,6 +496,41 @@ rule build_river_heat_potential: "../scripts/build_surface_water_heat_potentials/build_river_water_heat_potential.py" +rule build_lake_heat_potential: + params: + drop_leap_day=config_provider("enable", "drop_leap_day"), + snapshots=config_provider("snapshots"), + dh_area_buffer=config_provider( + "sector", "district_heating", "dh_areas", "buffer" + ), + enable_heat_source_maps=config_provider("plotting", "enable_heat_source_maps"), + input: + unpack(input_hera_data), + lake_data=rules.retrieve_lake_data_hera.output["lake_data"], + regions_onshore=resources("regions_onshore_base_s_{clusters}.geojson"), + dh_areas=resources("dh_areas_base_s_{clusters}.geojson"), + output: + heat_source_power=resources( + "heat_source_power_lake_water_base_s_{clusters}.csv" + ), + heat_source_temperature=resources("temp_lake_water_base_s_{clusters}.nc"), + heat_source_temperature_temporal_aggregate=resources( + "temp_lake_water_base_s_{clusters}_temporal_aggregate.nc" + ), + heat_source_energy_temporal_aggregate=resources( + "heat_source_energy_lake_water_base_s_{clusters}_temporal_aggregate.nc" + ), + resources: + mem_mb=20000, + log: + logs("build_lake_water_heat_potential_base_s_{clusters}.log"), + benchmark: + benchmarks("build_lake_water_heat_potential_base_s_{clusters}") + threads: 1 + script: + "../scripts/build_surface_water_heat_potentials/build_lake_water_heat_potential.py" + + def input_heat_source_temperature( w, replace_names: dict[str, str] = { diff --git a/rules/retrieve.smk b/rules/retrieve.smk index c1528660e3..a530e0979e 100755 --- a/rules/retrieve.smk +++ b/rules/retrieve.smk @@ -1158,3 +1158,33 @@ if (MOBILITY_PROFILES_DATASET := dataset_version("mobility_profiles"))["source"] run: copy2(input["kfz"], output["kfz"]) copy2(input["pkw"], output["pkw"]) + + +if (LAKE_DATA_DATASET := dataset_version("lake_data"))["source"] in [ + "primary", + "archive", +]: + + rule retrieve_lake_data: + input: + zip_file=storage(LAKE_DATA_DATASET["url"]), + output: + zip_file=f"{LAKE_DATA_DATASET['folder']}/HydroLAKES_polys_v10.gdb.zip", + lake_shapes=expand( + f"{LAKE_DATA_DATASET['folder']}/HydroLAKES_polys_v10.gdb.{{ext}}", + ext=[ + "shp", + "shx", + "dbf", + "cpg", + "prj", + "sbn", + "sbx", + ], + ), + run: + copy2(input["zip_file"], output["zip_file"]) + unpack_archive( + output["zip_file"], + LAKE_DATA_DATASET["folder"], + ) diff --git a/scripts/build_surface_water_heat_potentials/approximators/lake_water_heat_approximator.py b/scripts/build_surface_water_heat_potentials/approximators/lake_water_heat_approximator.py new file mode 100644 index 0000000000..839ef08696 --- /dev/null +++ b/scripts/build_surface_water_heat_potentials/approximators/lake_water_heat_approximator.py @@ -0,0 +1,114 @@ +# SPDX-FileCopyrightText: Contributors to PyPSA-Eur +# +# SPDX-License-Identifier: MIT +from functools import cached_property +from typing import Union + +import geopandas as gpd +import shapely +import xarray as xr + +from scripts.build_surface_water_heat_potentials.approximators.surface_water_heat_approximator import ( + RiverWaterHeatApproximator, + SurfaceWaterHeatApproximator, +) + + +class LakeWaterHeatApproximator(SurfaceWaterHeatApproximator): + """ + River water heat approximator for district heating systems. + + Parameters are mostly based on expert input and Triebs 2023: "Untersuchung der zukünftigen Fernwärmeversorgung unter Unsicherheit bei Berücksichtigung technischer, ökonomischer und ökologischer Randbedingungen". # codespell:ignore unter + + Min_distance of 25km is based on Jung et al.: "Estimation of + Temperature Recovery Distance and the Influence of Heat Pump Discharge on + Fluvial Ecosystems". + """ + + def __init__( + self, + ambient_temperature: xr.DataArray, + region: Union[shapely.geometry.polygon.Polygon, gpd.GeoSeries], + lake_shapes: gpd.GeoDataFrame, + delta_t_max: float = 1, + min_outlet_temperature: float = 1, + ) -> None: + self.ambient_temperature = ambient_temperature + self.lake_shapes = lake_shapes + + water_temperature = self._approximate_lake_temperature( + ambient_temperature=ambient_temperature + ) + water_temperature = water_temperature.rio.write_crs( + f"EPSG:{ambient_temperature.rio.crs.to_epsg()}" + ) + + super().__init__( + volume_flow=self._volume_flow_in_region, + water_temperature=self._round_coordinates(water_temperature), + region=region, + max_relative_volume_flow=1.0, + delta_t_max=delta_t_max, + min_outlet_temperature=min_outlet_temperature, + ) + + @staticmethod + def _round_coordinates(**kwargs) -> xr.DataArray: + return RiverWaterHeatApproximator._round_coordinates(**kwargs) + + @staticmethod + def _approximate_lake_temperature(**kwargs) -> xr.DataArray: + return RiverWaterHeatApproximator._approximate_river_temperature(**kwargs) + + @cached_property + def _volume_flow(self): + raise NotImplementedError() + + @cached_property + def _volume_flow_in_region(self): + lake_volume_km3 = self._lake_parts_in_regions.groupby("name")["Vol_total"].sum() + lake_volume_m3 = lake_volume_km3 * 1e9 + hourly_volume_m3 = lake_volume_m3 / 8760 + return hourly_volume_m3 + + @cached_property + def _lake_parts_in_region(self) -> gpd.GeoDataFrame: + """Get lake parts within the defined region.""" + lake_parts = gpd.overlay( + self.lake_shapes, + gpd.GeoDataFrame(geometry=[self.region], crs=self.lake_shapes.crs), + how="intersection", + ) + return lake_parts + + def _air_temperature_in_lakes( + self, eligible_lake_parts: gpd.GeoDataFrame + ) -> xr.DataArray: + """Get air temperature averaged over lake parts.""" + + return self.ambient_temperature.rio.clip( + eligible_lake_parts.geometry, eligible_lake_parts.crs, drop=True + ) + + @cached_property + def _water_temperature_in_region_raster(self) -> gpd.GeoDataFrame: + air_temperature_in_lakes = self._air_temperature_in_lakes( + self._lake_parts_in_region + ) + return self._approximate_lake_temperature(air_temperature_in_lakes) + + @cached_property + def _water_temperature_in_region(self) -> xr.DataArray: + return self._water_temperature_in_region_raster.mean(dim=("lat", "lon")) + + @cached_property + def _power_sum_spatial(self) -> None: + raise NotImplementedError("Lake power is not spatially resolved.") + + @cached_property + def _power_sum_temporal(self) -> None: + raise NotImplementedError("Lake power is not spatially resolved.") + + @cached_property + def _power_in_region(self): + raise NotImplementedError("Lake power is not spatially resolved.") diff --git a/scripts/build_surface_water_heat_potentials/approximators/river_water_heat_approximator.py b/scripts/build_surface_water_heat_potentials/approximators/river_water_heat_approximator.py index 5a9caa9424..2ba6bc5135 100644 --- a/scripts/build_surface_water_heat_potentials/approximators/river_water_heat_approximator.py +++ b/scripts/build_surface_water_heat_potentials/approximators/river_water_heat_approximator.py @@ -52,8 +52,9 @@ def __init__( min_distance_meters=min_distance_meters, ) + @staticmethod def _round_coordinates( - self, da: xr.DataArray, decimal_precision: int = 4 + da: xr.DataArray, decimal_precision: int = 4 ) -> xr.DataArray: """ Round the coordinates of the HERA dataset to the defined precision. diff --git a/scripts/build_surface_water_heat_potentials/build_lake_water_heat_potential.py b/scripts/build_surface_water_heat_potentials/build_lake_water_heat_potential.py new file mode 100644 index 0000000000..3bff71365b --- /dev/null +++ b/scripts/build_surface_water_heat_potentials/build_lake_water_heat_potential.py @@ -0,0 +1,444 @@ +# SPDX-FileCopyrightText: Contributors to PyPSA-Eur +# +# SPDX-License-Identifier: MIT +""" +Calculate river water heat potential for district heating systems. + +This script computes the thermal potential of rivers as a heat source for district +heating applications. It uses HERA river discharge and ambient temperature data to +estimate available heating power and average water temperatures across regions +intersected with district heating areas. + +The approximation accounts for temporal and spatial variations in river flow and temperature, +providing both spatial and temporal aggregates. Temporal aggregates are only used for plotting. + +Relevant Settings +----------------- + +.. code:: yaml + + sector: + district_heating: + dh_area_buffer: # Buffer around DH areas in meters to include nearby rivers + heat_source_cooling: # Exploitable temperature delta + snapshots: + start: + end: + enable: + drop_leap_day: + +Inputs +------ +- `data/hera_{year}/river_discharge_{year}.nc`: River discharge data from HERA +- `data/hera_{year}/ambient_temp_{year}.nc`: Ambient temperature data from HERA +- `resources//regions_onshore_base_s_{clusters}.geojson`: Onshore regions +- `resources//dh_areas_base_s_{clusters}.geojson`: District heating areas + +Outputs +------- +- `resources//heat_source_power_river_water_base_s_{clusters}.csv`: River heating power potentials by region +- `resources//temp_river_water_base_s_{clusters}.nc`: River water temperature profiles by region +- `resources//temp_river_water_base_s_{clusters}_temporal_aggregate.nc`: Temporal aggregated temperature data +- `resources//heat_source_energy_river_water_base_s_{clusters}_temporal_aggregate.nc`: Temporal aggregated energy data +""" + +import gc +import logging + +import dask +import geopandas as gpd +import numpy as np +import pandas as pd +import xarray as xr +from _helpers import ( + configure_logging, + get_snapshots, + set_scenario_config, + update_config_from_wildcards, +) +from approximators.river_water_heat_approximator import RiverWaterHeatApproximator + +logger = logging.getLogger(__name__) + +MEMORY_SAFETY_FACTOR = 0.7 # Use 70% of available memory for Dask arrays + + +def load_hera_data( + hera_inputs: dict, + snapshots: pd.DatetimeIndex, + minx: float, + miny: float, + maxx: float, + maxy: float, +) -> dict: + """ + Load and concatenate HERA data files for multiple years with spatial clipping. + + Parameters + ---------- + hera_inputs : dict + Dictionary with year-specific HERA file paths from input_hera_data(). + Expected keys: hera_river_discharge_{year}, hera_ambient_temperature_{year} + snapshots : pd.DatetimeIndex + Target snapshots to select from the combined data + minx, miny, maxx, maxy : float + Bounding box coordinates for spatial clipping + + Returns + ------- + dict + Dictionary with processed xarray DataArrays for river_discharge and ambient_temperature + """ + temp_files = [ + v for k, v in hera_inputs.items() if k.startswith("hera_ambient_temperature_") + ] + + # Determine time range from snapshots with buffer to ensure HERA coverage + # HERA data is on 6h intervals, so we need to extend the range to capture + # HERA timestamps that bracket our actual snapshots + buffer = pd.Timedelta(hours=12) # Buffer to ensure we capture surrounding HERA data + start_time = snapshots.min() - buffer + end_time = snapshots.max() + buffer + + # Load and concatenate ambient temperature files using open_mfdataset + ambient_temperature = xr.open_mfdataset( + temp_files, + chunks={"time": -1, "lat": 990, "lon": 1510}, + concat_dim="time", + combine="nested", + )["ta6"] + + # Select time range that covers our snapshots (using native HERA resolution) + ambient_temperature = ambient_temperature.sel(time=slice(start_time, end_time)) + + # Process ambient temperature data + return ( + ambient_temperature.rename({"lat": "latitude", "lon": "longitude"}) + .rio.write_crs("EPSG:4326") + .rio.clip_box(minx, miny, maxx, maxy) + .rio.reproject("EPSG:3035") + ) + + +def _create_empty_datasets( + snapshots: pd.DatetimeIndex, center_lon: float, center_lat: float +) -> tuple[xr.Dataset, xr.Dataset]: + """ + Create empty spatial and temporal aggregate datasets for regions without DH areas. + + When a region has no intersection with district heating areas, we still need + to provide valid datasets with zero values to maintain consistent data structure + across all regions. This prevents errors in downstream processing. + + Parameters + ---------- + snapshots : pd.DatetimeIndex + Time snapshots for the spatial aggregate + center_lon : float + Longitude of region center (for fallback coordinate) + center_lat : float + Latitude of region center (for fallback coordinate) + + Returns + ------- + tuple[xr.Dataset, xr.Dataset] + Tuple of (spatial_aggregate, temporal_aggregate) datasets with zero values + """ + # Create spatial aggregate with time-series of zeros + # This represents no available river heat power/temperature over time + spatial_aggregate = xr.Dataset( + data_vars={ + "total_power": xr.DataArray( + np.zeros(len(snapshots)), # Zero power for all timestamps + dims=["time"], + coords={"time": snapshots}, + ), + "average_temperature": xr.DataArray( + np.zeros(len(snapshots)), # Zero temperature for all timestamps + dims=["time"], + coords={"time": snapshots}, + ), + } + ) + + # Create temporal aggregate with single spatial point of zeros + # This represents no available river heat energy/temperature spatially + temporal_aggregate = xr.Dataset( + data_vars={ + "total_energy": xr.DataArray( + [[0.0]], # Zero energy at region center + dims=["longitude", "latitude"], + coords={"longitude": [center_lon], "latitude": [center_lat]}, + ), + "average_temperature": xr.DataArray( + [[0.0]], # Zero temperature at region center + dims=["longitude", "latitude"], + coords={"longitude": [center_lon], "latitude": [center_lat]}, + ), + } + ) + + return spatial_aggregate, temporal_aggregate + + +def get_regional_result( + hera_inputs: dict, + region: gpd.GeoSeries, + dh_areas: gpd.GeoDataFrame, + snapshots: pd.DatetimeIndex, + enable_heat_source_maps: bool = False, +) -> dict[str, xr.Dataset]: + """ + Calculate river water heat potential for a given region and district heating areas. + + Parameters + ---------- + hera_inputs : dict + Dictionary containing HERA input file paths with year-specific keys. + region : geopandas.GeoSeries + Geographical region for which to compute the heat potential. + dh_areas : geopandas.GeoDataFrame + District heating areas to intersect with the region. + snapshots : pd.DatetimeIndex + Time snapshots, used for loading data and for regions without dh_areas + enable_heat_source_maps : bool, default False + Whether to enable heat source mapping. + + Returns + ------- + dict + Dictionary with keys 'spatial aggregate' and 'temporal aggregate'. + 'spatial aggregate' contains total power and average temperature. + 'temporal aggregate' contains time series for energy and temperature (for analysis/plotting). + """ + # Store original region for fallback centroid calculation + original_region = region.copy() + + # Intersect region with district heating areas + intersected_geometry = gpd.overlay( + region.to_frame(), + dh_areas, + how="intersection", + ).union_all() + + region.geometry = intersected_geometry + + # Handle empty geometry case (no intersection with DH areas) + # This occurs when a region has no district heating area + # Note: the region could still have district heating, as central heat demand is computed elsewhere + if region.geometry.is_empty.any(): + # Get the center of the original region (before intersection) + # We use the original region to get a meaningful coordinate for the empty datasets + # Project to EPSG:3035 for accurate centroid calculation, then back to EPSG:4326 + region_center = ( + original_region.to_crs("EPSG:3035").centroid.to_crs("EPSG:4326").iloc[0] + ) + center_lon = region_center.x + center_lat = region_center.y + + # Return zero-filled datasets with proper structure + spatial_aggregate, temporal_aggregate = _create_empty_datasets( + snapshots, center_lon, center_lat + ) + + return { + "spatial aggregate": spatial_aggregate, + "temporal aggregate": temporal_aggregate, + } + + # Process region with valid DH area intersection + # Get bounding box for efficient data clipping + minx, miny, maxx, maxy = region.total_bounds + + # Data processing strategy: + # 1. Load HERA discharge and temperature data + # 2. Clip to region bounds for efficiency + # 3. Reproject to EPSG:3035 for accurate spatial calculations + # 4. Feed to approximator for heat potential calculation + + # Load and concatenate HERA data for all required years with preprocessing + hera_data = load_hera_data(hera_inputs, snapshots, minx, miny, maxx, maxy) + ambient_temperature = hera_data["ambient_temperature"] + + # Reproject region to match data CRS for spatial calculations + region = region.to_crs("EPSG:3035") + + lake_water_heat_approximator = RiverWaterHeatApproximator( + volume_flow=None, # River discharge (volume flow) + ambient_temperature=ambient_temperature, # Air temperature + region=region, # Geographic region of interest + ) + + # Calculate spatial aggregate (time series data for the entire region) + # Contains total_power [MW] and average_temperature [°C] over time + spatial_aggregate = lake_water_heat_approximator.get_spatial_aggregate() + + # Calculate temporal aggregate only if needed (spatial distribution data for plotting, no time dimension) + if enable_heat_source_maps: + temporal_aggregate = ( + lake_water_heat_approximator.get_temporal_aggregate() + .rio.reproject("EPSG:4326") # Convert back to WGS84 for output consistency + .rename({"x": "longitude", "y": "latitude"}) # Standardize coordinate names + ) + temporal_aggregate = temporal_aggregate.compute() + else: + temporal_aggregate = None + + # Compute results immediately to free Dask arrays and enable garbage collection + spatial_aggregate = spatial_aggregate.compute() + + # Explicitly delete approximator and intermediate arrays + del lake_water_heat_approximator + del ambient_temperature + gc.collect() + + result = { + "spatial aggregate": spatial_aggregate, + } + + # Only include temporal aggregate if computed + if temporal_aggregate is not None: + result["temporal aggregate"] = temporal_aggregate + if enable_heat_source_maps: + del temporal_aggregate + + return result + + +def set_dask_chunk_size( + n_threads: int, # Number of threads per worker, + memory_mb: int, # Memory per worker in MB + memory_safety_factor=MEMORY_SAFETY_FACTOR, + n_datasets: int = 2, # ambient temperature and river discharge datasets + operation_multiplier: int = 3, # Multiplier for operation overhead +) -> None: + """ + Set the Dask chunk size based on available memory and number of threads. + This function calculates the chunk size for Dask arrays to optimize memory usage + """ + + chunk_size = ( + memory_mb * memory_safety_factor / n_threads / n_datasets / operation_multiplier + ) + dask.config.set({"array.chunk-size": f"{chunk_size}MB"}) + + +if __name__ == "__main__": + if "snakemake" not in globals(): + from _helpers import mock_snakemake + + snakemake = mock_snakemake( + "build_river_water_heat_potential", + clusters="39", + opts="", + ll="vopt", + sector_opts="", + planning_horizons=2050, + ) + + # Configure logging and scenario + configure_logging(snakemake) + set_scenario_config(snakemake) + update_config_from_wildcards(snakemake.config, snakemake.wildcards) + + # Get simulation snapshots + snapshots: pd.DatetimeIndex = get_snapshots( + snakemake.params.snapshots, snakemake.params.drop_leap_day + ) + + # Load regions and district heating areas + regions_onshore = gpd.read_file(snakemake.input["regions_onshore"]) + regions_onshore.set_index("name", inplace=True) + regions_onshore = regions_onshore.to_crs("EPSG:4326") + + dh_areas = gpd.read_file(snakemake.input["dh_areas"]).to_crs("EPSG:3035") + # Buffer district heating areas by specified amount + dh_areas["geometry"] = dh_areas.geometry.buffer(snakemake.params.dh_area_buffer) + dh_areas = dh_areas.to_crs("EPSG:4326") + + lake_data = gpd.read_file(snakemake.input["lake_data"]).to_crs("EPSG:4326") + + # Configure Dask for multi-threading within operations (no distributed cluster) + dask.config.set(scheduler="threads") # Use threaded scheduler + dask.config.set(num_workers=snakemake.threads) # Use specified number of threads + + set_dask_chunk_size( + n_threads=snakemake.threads, memory_mb=snakemake.resources.mem_mb + ) + + # Process regions sequentially but with multi-threaded Dask operations + results = [] + for i, region_name in enumerate(regions_onshore.index, 1): + # Extract region geometry and create a copy to avoid modification conflicts + region = gpd.GeoSeries(regions_onshore.loc[region_name].copy(deep=True)) + + # Process region with multi-threaded Dask operations + result = get_regional_result( + hera_inputs=dict(snakemake.input), + region=region, + dh_areas=dh_areas, + snapshots=snapshots, + enable_heat_source_maps=snakemake.params.enable_heat_source_maps, + ) + results.append(result) + + # Explicit cleanup to free memory between regions + del result, region + gc.collect() + + # Build DataFrame of total power for each region + # Regions as columns and time as rows + power = pd.DataFrame( + { + region_name: res["spatial aggregate"]["total_power"].to_pandas() + for region_name, res in zip(regions_onshore.index, results) + } + ).dropna() + + power = power.reindex( + snapshots, method="nearest" + ) # Use "nearest" method to handle any minor timestamp differences due to floating point precision + + # Save power potentials in MW + power.to_csv(snakemake.output.heat_source_power) + + # Concatenate average temperature for all regions into single dataset + temperature = ( + xr.concat( + [res["spatial aggregate"]["average_temperature"] for res in results], + dim="name", + ) + .assign_coords(name=regions_onshore.index) + .dropna(dim="time") + ) + + # Align temperature data to snapshots, use nearest to handle any minor decimal differences + temperature = temperature.sel(time=snapshots, method="nearest").assign_coords( + time=snapshots + ) + + # Save temperature profiles as NetCDF for heat pump COP calculations + # Units: °C (degrees Celsius) + temperature.to_netcdf(snakemake.output.heat_source_temperature) + + # Save temporal aggregate results for analysis and visualization (if enabled) + if snakemake.params.enable_heat_source_maps: + # Energy temporal aggregate: spatial distribution of available energy + # Units: MWh (megawatt-hours) - total energy potential per location + xr.concat( + [res["temporal aggregate"]["total_energy"] for res in results], + dim=regions_onshore.index, + ).to_netcdf(snakemake.output.heat_source_energy_temporal_aggregate) + + # Temperature temporal aggregate: spatial distribution of temperatures + # Units: °C (degrees Celsius) - average temperature per location + xr.concat( + [res["temporal aggregate"]["average_temperature"] for res in results], + dim=regions_onshore.index, + ).to_netcdf(snakemake.output.heat_source_temperature_temporal_aggregate) + + else: + # Create empty placeholder files to satisfy Snakemake outputs + empty_ds = xr.Dataset() + empty_ds.to_netcdf(snakemake.output.heat_source_energy_temporal_aggregate) + empty_ds.to_netcdf(snakemake.output.heat_source_temperature_temporal_aggregate) From 24526ee8dae27e52ff2b9113b661b26757bdb26a Mon Sep 17 00:00:00 2001 From: Amos Schledorn Date: Mon, 22 Dec 2025 18:49:58 +0100 Subject: [PATCH 02/22] feat: update workflow --- config/config.default.yaml | 3 + data/versions.csv | 4 +- doc/data_inventory.csv | 3 +- rules/build_sector.smk | 2 +- rules/retrieve.smk | 13 +-- .../lake_water_heat_approximator.py | 85 +++++++++++++------ .../build_lake_water_heat_potential.py | 27 +++++- scripts/definitions/heat_source.py | 7 +- scripts/definitions/heat_system.py | 1 + 9 files changed, 98 insertions(+), 47 deletions(-) diff --git a/config/config.default.yaml b/config/config.default.yaml index 0ef3084c16..398445572e 100644 --- a/config/config.default.yaml +++ b/config/config.default.yaml @@ -1149,6 +1149,9 @@ data: jrc_ardeco: source: archive version: latest + lake_data: + source: archive + version: latest secrets: corine: '' #Add API key here if primary source is used for retrieving corine dataset after registering in https://land.copernicus.eu/user/login diff --git a/data/versions.csv b/data/versions.csv index a72c678ade..4d4f7ffc2a 100644 --- a/data/versions.csv +++ b/data/versions.csv @@ -99,5 +99,5 @@ "geothermal_heat_utilisation_potentials","archive","341.5","['latest', 'supported']","https://zenodo.org/records/17207640/files/Results_DH_Matching_Cluster.xlsx","" "jrc_ardeco","primary","2021","['latest', 'supported']","https://territorial.ec.europa.eu/ardeco-api-v2/rest/export/","" "jrc_ardeco","archive","2021","['latest', 'supported']","https://zenodo.org/records/17249457/files","" -"lake_data","primary","v10","['latest', 'supported']",""https://data.hydrosheds.org/file/hydrolakes/HydroLAKES_polys_v10.gdb.zip"","" -"lake_data","archive","v10","['latest', 'supported']",""https://data.hydrosheds.org/file/hydrolakes/HydroLAKES_polys_v10.gdb.zip"","" \ No newline at end of file +"lake_data","primary","v10","['latest', 'supported']","https://data.hydrosheds.org/file/hydrolakes/HydroLAKES_polys_v10.gdb.zip","" +"lake_data","archive","v10","['latest', 'supported']","https://data.hydrosheds.org/file/hydrolakes/HydroLAKES_polys_v10.gdb.zip","" \ No newline at end of file diff --git a/doc/data_inventory.csv b/doc/data_inventory.csv index 133ce39d6e..b3ffd818dd 100644 --- a/doc/data_inventory.csv +++ b/doc/data_inventory.csv @@ -41,4 +41,5 @@ "mobility_profiles","German Vehicle Activity Profiles","Vehicle activity profiles for different vehicle types and road types in Germany, based on monitoring data from the Federal Highway Research Institute (BASt). These profiles provide insights into travel behavior and patterns, which can be used for transport modeling and analysis.","Federal Highway Research Institute (BASt)","https://www.bast.de/DE/Themen/Digitales/HF_1/Massnahmen/verkehrszaehlung/Stundenwerte.html?nn=414410","CC-BY-4.0" "dh_areas","","Shapes of district heating areas","ISI Fraunhofer-Institut für System- und Innovationsforschung","https://fordatis.fraunhofer.de/handle/fordatis/341.5","CC-BY-4.0" "geothermal_heat_utilisation_potentials","","Potentials for Geothermal heat utilisation","ISI Fraunhofer-Institut für System- und Innovationsforschung","https://fordatis.fraunhofer.de/handle/fordatis/341.5","CC-BY-4.0" -"jrc_ardeco","Annual Regional Database of the European Commission's Directorate General for Regional and Urban Policy","The database contains a set of long time-series variables and indicators for EU regions, as well as for regions in some EFTA and candidate countries, at various statistical scales (NUTS1, NUTS2, NUTS3, metro regions).","European Commission","https://territorial.ec.europa.eu/ardeco","similar to CC-BY" \ No newline at end of file +"jrc_ardeco","Annual Regional Database of the European Commission's Directorate General for Regional and Urban Policy","The database contains a set of long time-series variables and indicators for EU regions, as well as for regions in some EFTA and candidate countries, at various statistical scales (NUTS1, NUTS2, NUTS3, metro regions).","European Commission","https://territorial.ec.europa.eu/ardeco","similar to CC-BY" +"lake_data","HydroLakes database","Shapes, volumes and further data on global lakes","Messager, M.L., Lehner, B., Grill, G., Nedeva, I., Schmitt, O. (2016)","https://www.hydrosheds.org/products/hydrolakes","CC-BY-4.0" \ No newline at end of file diff --git a/rules/build_sector.smk b/rules/build_sector.smk index c519c121df..678461b725 100755 --- a/rules/build_sector.smk +++ b/rules/build_sector.smk @@ -506,7 +506,7 @@ rule build_lake_heat_potential: enable_heat_source_maps=config_provider("plotting", "enable_heat_source_maps"), input: unpack(input_hera_data), - lake_data=rules.retrieve_lake_data_hera.output["lake_data"], + lake_data=rules.retrieve_lake_data.output["lake_data"], regions_onshore=resources("regions_onshore_base_s_{clusters}.geojson"), dh_areas=resources("dh_areas_base_s_{clusters}.geojson"), output: diff --git a/rules/retrieve.smk b/rules/retrieve.smk index a530e0979e..09d4936332 100755 --- a/rules/retrieve.smk +++ b/rules/retrieve.smk @@ -1170,17 +1170,8 @@ if (LAKE_DATA_DATASET := dataset_version("lake_data"))["source"] in [ zip_file=storage(LAKE_DATA_DATASET["url"]), output: zip_file=f"{LAKE_DATA_DATASET['folder']}/HydroLAKES_polys_v10.gdb.zip", - lake_shapes=expand( - f"{LAKE_DATA_DATASET['folder']}/HydroLAKES_polys_v10.gdb.{{ext}}", - ext=[ - "shp", - "shx", - "dbf", - "cpg", - "prj", - "sbn", - "sbx", - ], + lake_data=directory( + f"{LAKE_DATA_DATASET['folder']}/HydroLAKES_polys_v10.gdb" ), run: copy2(input["zip_file"], output["zip_file"]) diff --git a/scripts/build_surface_water_heat_potentials/approximators/lake_water_heat_approximator.py b/scripts/build_surface_water_heat_potentials/approximators/lake_water_heat_approximator.py index 839ef08696..0e4906120d 100644 --- a/scripts/build_surface_water_heat_potentials/approximators/lake_water_heat_approximator.py +++ b/scripts/build_surface_water_heat_potentials/approximators/lake_water_heat_approximator.py @@ -8,8 +8,10 @@ import shapely import xarray as xr -from scripts.build_surface_water_heat_potentials.approximators.surface_water_heat_approximator import ( +from scripts.build_surface_water_heat_potentials.approximators.river_water_heat_approximator import ( RiverWaterHeatApproximator, +) +from scripts.build_surface_water_heat_potentials.approximators.surface_water_heat_approximator import ( SurfaceWaterHeatApproximator, ) @@ -35,6 +37,7 @@ def __init__( ) -> None: self.ambient_temperature = ambient_temperature self.lake_shapes = lake_shapes + self.region = region # Set early so _lake_parts_in_region can access it water_temperature = self._approximate_lake_temperature( ambient_temperature=ambient_temperature @@ -53,30 +56,69 @@ def __init__( ) @staticmethod - def _round_coordinates(**kwargs) -> xr.DataArray: - return RiverWaterHeatApproximator._round_coordinates(**kwargs) + def _round_coordinates( + da: xr.DataArray, decimal_precision: int = 4 + ) -> xr.DataArray: + return RiverWaterHeatApproximator._round_coordinates(da, decimal_precision) @staticmethod def _approximate_lake_temperature(**kwargs) -> xr.DataArray: return RiverWaterHeatApproximator._approximate_river_temperature(**kwargs) @cached_property - def _volume_flow(self): - raise NotImplementedError() + def _total_volume_flow_m3_per_s(self) -> float: + """ + Calculate total volume flow in m³/s from lake volume. + + Vol_total in HydroLAKES is in MCM (million cubic meters = 1e6 m³). + """ + lake_volume_mcm = self._lake_parts_in_region["Vol_total"].sum() + lake_volume_m3 = lake_volume_mcm * 1e6 # MCM to m³ + # Convert yearly turnover to per-second flow (m³/s) + volume_flow_m3_per_s = lake_volume_m3 / (8760 * 3600) + return volume_flow_m3_per_s @cached_property - def _volume_flow_in_region(self): - lake_volume_km3 = self._lake_parts_in_regions.groupby("name")["Vol_total"].sum() - lake_volume_m3 = lake_volume_km3 * 1e9 - hourly_volume_m3 = lake_volume_m3 / 8760 - return hourly_volume_m3 + def _volume_flow_in_region(self) -> xr.DataArray: + """ + Create a volume flow raster distributed uniformly over lake pixels. + + The total flow is divided by the number of valid pixels so that + sum(flow_raster) = total_flow. + """ + total_flow = self._total_volume_flow_m3_per_s + + # Use water temperature raster as template for the shape (before averaging) + template = self._water_temperature_in_region_raster + + # Count valid (non-NaN) pixels in the raster + # Use a 2D slice to count pixels (same count for all time steps) + template_2d = template.isel(time=0) if "time" in template.dims else template + n_valid_pixels = template_2d.notnull().sum().values + + if n_valid_pixels == 0: + return xr.full_like(template, 0.0) + + # Distribute flow uniformly so sum = total_flow + flow_per_pixel = total_flow / n_valid_pixels + + # Create raster with flow_per_pixel where template is valid, 0 elsewhere + flow_raster = xr.where(template.notnull(), flow_per_pixel, 0.0) + return flow_raster @cached_property def _lake_parts_in_region(self) -> gpd.GeoDataFrame: """Get lake parts within the defined region.""" + # Handle both GeoSeries and single geometry + if isinstance(self.region, gpd.GeoSeries): + region_gdf = gpd.GeoDataFrame(geometry=self.region, crs=self.region.crs) + else: + region_gdf = gpd.GeoDataFrame( + geometry=[self.region], crs=self.lake_shapes.crs + ) lake_parts = gpd.overlay( self.lake_shapes, - gpd.GeoDataFrame(geometry=[self.region], crs=self.lake_shapes.crs), + region_gdf, how="intersection", ) return lake_parts @@ -91,24 +133,15 @@ def _air_temperature_in_lakes( ) @cached_property - def _water_temperature_in_region_raster(self) -> gpd.GeoDataFrame: + def _water_temperature_in_region_raster(self) -> xr.DataArray: air_temperature_in_lakes = self._air_temperature_in_lakes( self._lake_parts_in_region ) - return self._approximate_lake_temperature(air_temperature_in_lakes) + return self._approximate_lake_temperature( + ambient_temperature=air_temperature_in_lakes + ) @cached_property def _water_temperature_in_region(self) -> xr.DataArray: - return self._water_temperature_in_region_raster.mean(dim=("lat", "lon")) - - @cached_property - def _power_sum_spatial(self) -> None: - raise NotImplementedError("Lake power is not spatially resolved.") - - @cached_property - def _power_sum_temporal(self) -> None: - raise NotImplementedError("Lake power is not spatially resolved.") - - @cached_property - def _power_in_region(self): - raise NotImplementedError("Lake power is not spatially resolved.") + # Use x, y dims (EPSG:3035 projected) instead of lat, lon + return self._water_temperature_in_region_raster.mean(dim=("x", "y")) diff --git a/scripts/build_surface_water_heat_potentials/build_lake_water_heat_potential.py b/scripts/build_surface_water_heat_potentials/build_lake_water_heat_potential.py index 3bff71365b..653fdab7e8 100644 --- a/scripts/build_surface_water_heat_potentials/build_lake_water_heat_potential.py +++ b/scripts/build_surface_water_heat_potentials/build_lake_water_heat_potential.py @@ -56,7 +56,7 @@ set_scenario_config, update_config_from_wildcards, ) -from approximators.river_water_heat_approximator import RiverWaterHeatApproximator +from approximators.lake_water_heat_approximator import LakeWaterHeatApproximator logger = logging.getLogger(__name__) @@ -112,13 +112,15 @@ def load_hera_data( ambient_temperature = ambient_temperature.sel(time=slice(start_time, end_time)) # Process ambient temperature data - return ( + ambient_temperature = ( ambient_temperature.rename({"lat": "latitude", "lon": "longitude"}) .rio.write_crs("EPSG:4326") .rio.clip_box(minx, miny, maxx, maxy) .rio.reproject("EPSG:3035") ) + return {"ambient_temperature": ambient_temperature} + def _create_empty_datasets( snapshots: pd.DatetimeIndex, center_lon: float, center_lat: float @@ -185,6 +187,7 @@ def get_regional_result( hera_inputs: dict, region: gpd.GeoSeries, dh_areas: gpd.GeoDataFrame, + lake_shapes: gpd.GeoDataFrame, snapshots: pd.DatetimeIndex, enable_heat_source_maps: bool = False, ) -> dict[str, xr.Dataset]: @@ -262,11 +265,12 @@ def get_regional_result( # Reproject region to match data CRS for spatial calculations region = region.to_crs("EPSG:3035") + lake_shapes = lake_shapes.to_crs("EPSG:3035") - lake_water_heat_approximator = RiverWaterHeatApproximator( - volume_flow=None, # River discharge (volume flow) + lake_water_heat_approximator = LakeWaterHeatApproximator( ambient_temperature=ambient_temperature, # Air temperature region=region, # Geographic region of interest + lake_shapes=lake_shapes, # Lake polygons with Vol_total ) # Calculate spatial aggregate (time series data for the entire region) @@ -377,6 +381,7 @@ def set_dask_chunk_size( hera_inputs=dict(snakemake.input), region=region, dh_areas=dh_areas, + lake_shapes=lake_data, snapshots=snapshots, enable_heat_source_maps=snakemake.params.enable_heat_source_maps, ) @@ -402,6 +407,13 @@ def set_dask_chunk_size( # Save power potentials in MW power.to_csv(snakemake.output.heat_source_power) + # Log computed power and temperature values + power_mean = power.mean() + logger.info("=== Lake heat potential summary ===") + logger.info("Mean power per region (MW):") + for region in power_mean.index: + logger.info(f" {region}: {power_mean[region]:.2f} MW") + # Concatenate average temperature for all regions into single dataset temperature = ( xr.concat( @@ -421,6 +433,13 @@ def set_dask_chunk_size( # Units: °C (degrees Celsius) temperature.to_netcdf(snakemake.output.heat_source_temperature) + # Log temperature values + temp_mean = temperature.mean(dim="time").to_pandas() + logger.info("Mean temperature per region (°C):") + for region in temp_mean.index: + logger.info(f" {region}: {temp_mean[region]:.2f} °C") + logger.info("=== END summary ===") + # Save temporal aggregate results for analysis and visualization (if enabled) if snakemake.params.enable_heat_source_maps: # Energy temporal aggregate: spatial distribution of available energy diff --git a/scripts/definitions/heat_source.py b/scripts/definitions/heat_source.py index 4534982571..7ac743c85f 100644 --- a/scripts/definitions/heat_source.py +++ b/scripts/definitions/heat_source.py @@ -42,7 +42,7 @@ class HeatSource(Enum): **Inexhaustible sources** (AIR, GROUND, SEA_WATER): Always available, no resource bus needed, heat pump draws from ambient. - **Limited sources requiring a bus** (GEOTHERMAL, RIVER_WATER, PTES): + **Limited sources requiring a bus** (GEOTHERMAL, RIVER_WATER, LAKE_WATER, PTES): Have spatial/temporal constraints, require resource tracking via buses. May support direct utilisation or preheating depending on temperature. @@ -58,6 +58,8 @@ class HeatSource(Enum): River water heat source with time-varying temperature. SEA_WATER : str Sea water heat source (treated as inexhaustible). + LAKE_WATER : str + Lake water heat source (treated as inexhaustible). AIR : str Ambient air heat source (inexhaustible). GROUND : str @@ -75,6 +77,7 @@ class HeatSource(Enum): GEOTHERMAL = "geothermal" RIVER_WATER = "river_water" SEA_WATER = "sea_water" + LAKE_WATER = "lake_water" AIR = "air" GROUND = "ground" PTES = "ptes" @@ -130,7 +133,7 @@ def temperature_from_config(self) -> bool: True for sources with config-defined temperatures (geothermal, PTX). False for sources with file-based time-series (river_water, ptes). """ - if self == HeatSource.RIVER_WATER: + if self in [HeatSource.RIVER_WATER, HeatSource.LAKE_WATER]: return False return self.source_type in [ HeatSourceType.SUPPLY_LIMITED, diff --git a/scripts/definitions/heat_system.py b/scripts/definitions/heat_system.py index 3b797927f2..b63d031527 100644 --- a/scripts/definitions/heat_system.py +++ b/scripts/definitions/heat_system.py @@ -234,6 +234,7 @@ def heat_pump_costs_name(self, heat_source: HeatSource | str) -> str: HeatSource.GEOTHERMAL, HeatSource.SEA_WATER, HeatSource.RIVER_WATER, + HeatSource.LAKE_WATER, HeatSource.ELECTROLYSIS_EXCESS, HeatSource.FISCHER_TROPSCH_EXCESS, HeatSource.SABATIER_EXCESS, From 9e9d3dc09ec0e7f1ffe66734f5cabff9865e4159 Mon Sep 17 00:00:00 2001 From: Amos Schledorn Date: Tue, 23 Dec 2025 12:28:20 +0100 Subject: [PATCH 03/22] remove obsolete comment --- .../approximators/surface_water_heat_approximator.py | 4 ---- 1 file changed, 4 deletions(-) diff --git a/scripts/build_surface_water_heat_potentials/approximators/surface_water_heat_approximator.py b/scripts/build_surface_water_heat_potentials/approximators/surface_water_heat_approximator.py index ed4bb4bc32..38acf5bb84 100644 --- a/scripts/build_surface_water_heat_potentials/approximators/surface_water_heat_approximator.py +++ b/scripts/build_surface_water_heat_potentials/approximators/surface_water_heat_approximator.py @@ -187,10 +187,6 @@ def _validate_and_reproject_input(self) -> None: f"{coord} coordinate '{coord}' not found in both datasets" ) - # For region geometry, we just check the type - # if not isinstance(self.region, shapely.geometry.multipolygon.MultiPolygon): - # raise ValueError(f"region_geometry must be a shapely MultiPolygon, got {type(self.region)}") - @cached_property def _volume_flow_in_region(self) -> xr.DataArray: """ From 26f3877c7434ecc19d1bf559f4ebe9ede934a463 Mon Sep 17 00:00:00 2001 From: Amos Schledorn Date: Tue, 23 Dec 2025 12:31:38 +0100 Subject: [PATCH 04/22] fix lake heat computations --- .../lake_water_heat_approximator.py | 90 ++++++++++--------- .../build_lake_water_heat_potential.py | 20 ++--- 2 files changed, 54 insertions(+), 56 deletions(-) diff --git a/scripts/build_surface_water_heat_potentials/approximators/lake_water_heat_approximator.py b/scripts/build_surface_water_heat_potentials/approximators/lake_water_heat_approximator.py index 0e4906120d..d7f9888875 100644 --- a/scripts/build_surface_water_heat_potentials/approximators/lake_water_heat_approximator.py +++ b/scripts/build_surface_water_heat_potentials/approximators/lake_water_heat_approximator.py @@ -33,10 +33,17 @@ def __init__( region: Union[shapely.geometry.polygon.Polygon, gpd.GeoSeries], lake_shapes: gpd.GeoDataFrame, delta_t_max: float = 1, - min_outlet_temperature: float = 1, + min_outlet_temperature: float = -100, ) -> None: self.ambient_temperature = ambient_temperature + self.lake_shapes = lake_shapes + lake_shapes.geometry = lake_shapes.geometry.map( + lambda g: type(g)([p.exterior for p in g.geoms]) + ).simplify( # drop holes + 0.001, preserve_topology=True + ) # simplify (~100 m) + self.region = region # Set early so _lake_parts_in_region can access it water_temperature = self._approximate_lake_temperature( @@ -48,63 +55,33 @@ def __init__( super().__init__( volume_flow=self._volume_flow_in_region, - water_temperature=self._round_coordinates(water_temperature), + water_temperature=water_temperature, region=region, max_relative_volume_flow=1.0, delta_t_max=delta_t_max, min_outlet_temperature=min_outlet_temperature, ) - @staticmethod - def _round_coordinates( - da: xr.DataArray, decimal_precision: int = 4 - ) -> xr.DataArray: - return RiverWaterHeatApproximator._round_coordinates(da, decimal_precision) + @property + def _scaling_factor(self) -> float: + return 1.0 # No scaling needed for lakes @staticmethod def _approximate_lake_temperature(**kwargs) -> xr.DataArray: return RiverWaterHeatApproximator._approximate_river_temperature(**kwargs) @cached_property - def _total_volume_flow_m3_per_s(self) -> float: + def _volume_flow_in_region(self) -> float: """ Calculate total volume flow in m³/s from lake volume. Vol_total in HydroLAKES is in MCM (million cubic meters = 1e6 m³). """ - lake_volume_mcm = self._lake_parts_in_region["Vol_total"].sum() - lake_volume_m3 = lake_volume_mcm * 1e6 # MCM to m³ - # Convert yearly turnover to per-second flow (m³/s) - volume_flow_m3_per_s = lake_volume_m3 / (8760 * 3600) - return volume_flow_m3_per_s - - @cached_property - def _volume_flow_in_region(self) -> xr.DataArray: - """ - Create a volume flow raster distributed uniformly over lake pixels. - - The total flow is divided by the number of valid pixels so that - sum(flow_raster) = total_flow. - """ - total_flow = self._total_volume_flow_m3_per_s - - # Use water temperature raster as template for the shape (before averaging) - template = self._water_temperature_in_region_raster - - # Count valid (non-NaN) pixels in the raster - # Use a 2D slice to count pixels (same count for all time steps) - template_2d = template.isel(time=0) if "time" in template.dims else template - n_valid_pixels = template_2d.notnull().sum().values - - if n_valid_pixels == 0: - return xr.full_like(template, 0.0) - - # Distribute flow uniformly so sum = total_flow - flow_per_pixel = total_flow / n_valid_pixels - - # Create raster with flow_per_pixel where template is valid, 0 elsewhere - flow_raster = xr.where(template.notnull(), flow_per_pixel, 0.0) - return flow_raster + lake_volume_km3 = self._lake_parts_in_region["Vol_total"].sum() + lake_volume_m3 = lake_volume_km3 * 1e6 # MCM to m³ + # Convert yearly turnover to per-hour flow (m³/h) + volume_flow_m3_per_h = lake_volume_m3 / 8760 + return volume_flow_m3_per_h @cached_property def _lake_parts_in_region(self) -> gpd.GeoDataFrame: @@ -143,5 +120,34 @@ def _water_temperature_in_region_raster(self) -> xr.DataArray: @cached_property def _water_temperature_in_region(self) -> xr.DataArray: - # Use x, y dims (EPSG:3035 projected) instead of lat, lon return self._water_temperature_in_region_raster.mean(dim=("x", "y")) + + @cached_property + def _power_sum_spatial(self) -> xr.DataArray: + """ + Cache the expensive spatial sum of power. + + Returns + ------- + xr.DataArray + Spatial sum of power over x and y dimensions + """ + return self._power_in_region + + def get_spatial_aggregate(self) -> xr.Dataset: + """ + Get the spatial aggregate of water temperature and power. + + Returns + ------- + xr.Dataset + Dataset containing total_power and average_temperature + """ + # Data is already spatially aggregated + # Combine into a single dataset + return xr.Dataset( + data_vars={ + "total_power": self._power_sum_spatial, + "average_temperature": self._water_temperature_in_region, + } + ) diff --git a/scripts/build_surface_water_heat_potentials/build_lake_water_heat_potential.py b/scripts/build_surface_water_heat_potentials/build_lake_water_heat_potential.py index 653fdab7e8..66f9dc2886 100644 --- a/scripts/build_surface_water_heat_potentials/build_lake_water_heat_potential.py +++ b/scripts/build_surface_water_heat_potentials/build_lake_water_heat_potential.py @@ -112,15 +112,13 @@ def load_hera_data( ambient_temperature = ambient_temperature.sel(time=slice(start_time, end_time)) # Process ambient temperature data - ambient_temperature = ( + return ( ambient_temperature.rename({"lat": "latitude", "lon": "longitude"}) .rio.write_crs("EPSG:4326") .rio.clip_box(minx, miny, maxx, maxy) .rio.reproject("EPSG:3035") ) - return {"ambient_temperature": ambient_temperature} - def _create_empty_datasets( snapshots: pd.DatetimeIndex, center_lon: float, center_lat: float @@ -253,15 +251,8 @@ def get_regional_result( # Get bounding box for efficient data clipping minx, miny, maxx, maxy = region.total_bounds - # Data processing strategy: - # 1. Load HERA discharge and temperature data - # 2. Clip to region bounds for efficiency - # 3. Reproject to EPSG:3035 for accurate spatial calculations - # 4. Feed to approximator for heat potential calculation - # Load and concatenate HERA data for all required years with preprocessing - hera_data = load_hera_data(hera_inputs, snapshots, minx, miny, maxx, maxy) - ambient_temperature = hera_data["ambient_temperature"] + ambient_temperature = load_hera_data(hera_inputs, snapshots, minx, miny, maxx, maxy) # Reproject region to match data CRS for spatial calculations region = region.to_crs("EPSG:3035") @@ -280,9 +271,10 @@ def get_regional_result( # Calculate temporal aggregate only if needed (spatial distribution data for plotting, no time dimension) if enable_heat_source_maps: temporal_aggregate = ( - lake_water_heat_approximator.get_temporal_aggregate() - .rio.reproject("EPSG:4326") # Convert back to WGS84 for output consistency - .rename({"x": "longitude", "y": "latitude"}) # Standardize coordinate names + lake_water_heat_approximator.get_temporal_aggregate().rio.reproject( + "EPSG:4326" + ) # Convert back to WGS84 for output consistency + # .rename({"x": "longitude", "y": "latitude"}) # Standardize coordinate names ) temporal_aggregate = temporal_aggregate.compute() else: From 3d8490233d2a37c6680b2ae82ca28b7479460a91 Mon Sep 17 00:00:00 2001 From: Amos Schledorn Date: Tue, 23 Dec 2025 14:26:09 +0100 Subject: [PATCH 05/22] update docs --- .../lake_water_heat_approximator.py | 176 +++++++++++++---- .../build_lake_water_heat_potential.py | 187 +++++++++--------- 2 files changed, 230 insertions(+), 133 deletions(-) diff --git a/scripts/build_surface_water_heat_potentials/approximators/lake_water_heat_approximator.py b/scripts/build_surface_water_heat_potentials/approximators/lake_water_heat_approximator.py index d7f9888875..9ce85eae7a 100644 --- a/scripts/build_surface_water_heat_potentials/approximators/lake_water_heat_approximator.py +++ b/scripts/build_surface_water_heat_potentials/approximators/lake_water_heat_approximator.py @@ -1,6 +1,8 @@ # SPDX-FileCopyrightText: Contributors to PyPSA-Eur # # SPDX-License-Identifier: MIT +"""Lake water heat approximator for district heating systems.""" + from functools import cached_property from typing import Union @@ -18,13 +20,34 @@ class LakeWaterHeatApproximator(SurfaceWaterHeatApproximator): """ - River water heat approximator for district heating systems. + Lake water heat approximator for district heating systems. + + Estimates the thermal potential of lakes as a heat source for district + heating applications. Uses lake volume data from HydroLAKES and ambient + temperature from HERA to compute available heating power. - Parameters are mostly based on expert input and Triebs 2023: "Untersuchung der zukünftigen Fernwärmeversorgung unter Unsicherheit bei Berücksichtigung technischer, ökonomischer und ökologischer Randbedingungen". # codespell:ignore unter + The water temperature approximation follows the methodology from + Triebs & Tsatsaronis 2022 (originally developed for rivers). - Min_distance of 25km is based on Jung et al.: "Estimation of - Temperature Recovery Distance and the Influence of Heat Pump Discharge on - Fluvial Ecosystems". + Parameters + ---------- + ambient_temperature : xr.DataArray + Ambient air temperature data (HERA), used to approximate water temperature. + region : Union[shapely.geometry.polygon.Polygon, gpd.GeoSeries] + Geographic region of interest. + lake_shapes : gpd.GeoDataFrame + Lake polygons from HydroLAKES with Vol_total attribute. + delta_t_max : float, optional + Maximum temperature difference for heat extraction, by default 1 K. + min_outlet_temperature : float, optional + Minimum outlet water temperature, by default 1 °C. + + References + ---------- + .. [1] Messager et al. (2016). Estimating the volume and age of water + stored in global lakes. Nature Communications, 13603. + .. [2] Triebs & Tsatsaronis (2022). Estimating the local renewable + potentials for the transformation of district heating systems. """ def __init__( @@ -33,18 +56,11 @@ def __init__( region: Union[shapely.geometry.polygon.Polygon, gpd.GeoSeries], lake_shapes: gpd.GeoDataFrame, delta_t_max: float = 1, - min_outlet_temperature: float = -100, + min_outlet_temperature: float = 1, ) -> None: self.ambient_temperature = ambient_temperature - - self.lake_shapes = lake_shapes - lake_shapes.geometry = lake_shapes.geometry.map( - lambda g: type(g)([p.exterior for p in g.geoms]) - ).simplify( # drop holes - 0.001, preserve_topology=True - ) # simplify (~100 m) - - self.region = region # Set early so _lake_parts_in_region can access it + self.lake_shapes = self._simplify_lake_geometries(lake_shapes) + self.region = region water_temperature = self._approximate_lake_temperature( ambient_temperature=ambient_temperature @@ -64,53 +80,135 @@ def __init__( @property def _scaling_factor(self) -> float: - return 1.0 # No scaling needed for lakes + """Scaling factor for power calculation (disabled for lakes).""" + return 1.0 + + @staticmethod + def _simplify_lake_geometries( + lake_shapes: gpd.GeoDataFrame, tolerance: float = 0.001 + ) -> gpd.GeoDataFrame: + """ + Simplify lake polygon geometries for faster spatial operations. + + Removes interior holes and simplifies exterior boundaries. + + Parameters + ---------- + lake_shapes : gpd.GeoDataFrame + Lake polygons from HydroLAKES. + tolerance : float, optional + Simplification tolerance in CRS units (~100m for EPSG:3035), + by default 0.001. + + Returns + ------- + gpd.GeoDataFrame + Lake shapes with simplified geometries. + """ + from shapely.geometry import MultiPolygon, Polygon + + def remove_holes(geom): + """Remove interior holes from polygon geometries.""" + if isinstance(geom, Polygon): + return Polygon(geom.exterior) + elif isinstance(geom, MultiPolygon): + return MultiPolygon([Polygon(p.exterior) for p in geom.geoms]) + return geom + + lake_shapes = lake_shapes.copy() + lake_shapes.geometry = lake_shapes.geometry.map(remove_holes).simplify( + tolerance, preserve_topology=True + ) + return lake_shapes @staticmethod def _approximate_lake_temperature(**kwargs) -> xr.DataArray: + """ + Approximate lake water temperature from ambient air temperature. + + Uses the river temperature approximation from Triebs & Tsatsaronis 2022, + which applies a moving average and empirical formula to derive water + temperature from ambient air temperature. + + Parameters + ---------- + **kwargs + Keyword arguments passed to + ``RiverWaterHeatApproximator._approximate_river_temperature``. + + Returns + ------- + xr.DataArray + Approximated lake water temperature. + """ return RiverWaterHeatApproximator._approximate_river_temperature(**kwargs) @cached_property def _volume_flow_in_region(self) -> float: """ - Calculate total volume flow in m³/s from lake volume. + Calculate volume flow rate from total lake volume. + Assumes annual turnover of lake volume, converted to hourly flow rate. Vol_total in HydroLAKES is in MCM (million cubic meters = 1e6 m³). + + Returns + ------- + float + Volume flow rate in m³/h. """ - lake_volume_km3 = self._lake_parts_in_region["Vol_total"].sum() - lake_volume_m3 = lake_volume_km3 * 1e6 # MCM to m³ - # Convert yearly turnover to per-hour flow (m³/h) + lake_volume_mcm = self._lake_parts_in_region["Vol_total"].sum() + lake_volume_m3 = lake_volume_mcm * 1e6 volume_flow_m3_per_h = lake_volume_m3 / 8760 return volume_flow_m3_per_h @cached_property def _lake_parts_in_region(self) -> gpd.GeoDataFrame: - """Get lake parts within the defined region.""" - # Handle both GeoSeries and single geometry + """ + Get lake polygons intersected with the region. + + Returns + ------- + gpd.GeoDataFrame + Lake parts that fall within the defined region. + """ if isinstance(self.region, gpd.GeoSeries): region_gdf = gpd.GeoDataFrame(geometry=self.region, crs=self.region.crs) else: region_gdf = gpd.GeoDataFrame( geometry=[self.region], crs=self.lake_shapes.crs ) - lake_parts = gpd.overlay( - self.lake_shapes, - region_gdf, - how="intersection", - ) - return lake_parts + return gpd.overlay(self.lake_shapes, region_gdf, how="intersection") def _air_temperature_in_lakes( self, eligible_lake_parts: gpd.GeoDataFrame ) -> xr.DataArray: - """Get air temperature averaged over lake parts.""" + """ + Clip ambient temperature raster to lake areas. + + Parameters + ---------- + eligible_lake_parts : gpd.GeoDataFrame + Lake polygons to clip temperature data to. + Returns + ------- + xr.DataArray + Ambient temperature clipped to lake areas. + """ return self.ambient_temperature.rio.clip( eligible_lake_parts.geometry, eligible_lake_parts.crs, drop=True ) @cached_property def _water_temperature_in_region_raster(self) -> xr.DataArray: + """ + Compute lake water temperature raster for the region. + + Returns + ------- + xr.DataArray + Water temperature raster with spatial (x, y) and time dimensions. + """ air_temperature_in_lakes = self._air_temperature_in_lakes( self._lake_parts_in_region ) @@ -120,31 +218,39 @@ def _water_temperature_in_region_raster(self) -> xr.DataArray: @cached_property def _water_temperature_in_region(self) -> xr.DataArray: + """ + Compute spatially-averaged lake water temperature. + + Returns + ------- + xr.DataArray + Mean water temperature over time (spatial dimensions averaged). + """ return self._water_temperature_in_region_raster.mean(dim=("x", "y")) @cached_property def _power_sum_spatial(self) -> xr.DataArray: """ - Cache the expensive spatial sum of power. + Get total power (already spatially aggregated for lakes). Returns ------- xr.DataArray - Spatial sum of power over x and y dimensions + Total thermal power over time. """ return self._power_in_region def get_spatial_aggregate(self) -> xr.Dataset: """ - Get the spatial aggregate of water temperature and power. + Get spatially aggregated heat potential results. Returns ------- xr.Dataset - Dataset containing total_power and average_temperature + Dataset with variables: + - total_power : Total thermal power [MW] over time. + - average_temperature : Mean water temperature [°C] over time. """ - # Data is already spatially aggregated - # Combine into a single dataset return xr.Dataset( data_vars={ "total_power": self._power_sum_spatial, diff --git a/scripts/build_surface_water_heat_potentials/build_lake_water_heat_potential.py b/scripts/build_surface_water_heat_potentials/build_lake_water_heat_potential.py index 66f9dc2886..bc2ca1a8b2 100644 --- a/scripts/build_surface_water_heat_potentials/build_lake_water_heat_potential.py +++ b/scripts/build_surface_water_heat_potentials/build_lake_water_heat_potential.py @@ -2,15 +2,16 @@ # # SPDX-License-Identifier: MIT """ -Calculate river water heat potential for district heating systems. +Calculate lake water heat potential for district heating systems. -This script computes the thermal potential of rivers as a heat source for district -heating applications. It uses HERA river discharge and ambient temperature data to -estimate available heating power and average water temperatures across regions +This script computes the thermal potential of lakes as a heat source for district +heating applications. It uses HydroLAKES lake volume data and HERA ambient temperature +data to estimate available heating power and average water temperatures across regions intersected with district heating areas. -The approximation accounts for temporal and spatial variations in river flow and temperature, -providing both spatial and temporal aggregates. Temporal aggregates are only used for plotting. +The approximation accounts for temporal and spatial variations in lake volume and +temperature, providing both spatial and temporal aggregates. Temporal aggregates are +only used for plotting. Relevant Settings ----------------- @@ -19,8 +20,10 @@ sector: district_heating: - dh_area_buffer: # Buffer around DH areas in meters to include nearby rivers + dh_area_buffer: # Buffer around DH areas in meters to include nearby lakes heat_source_cooling: # Exploitable temperature delta + heat_sources: + - lake_water snapshots: start: end: @@ -29,21 +32,22 @@ Inputs ------ -- `data/hera_{year}/river_discharge_{year}.nc`: River discharge data from HERA -- `data/hera_{year}/ambient_temp_{year}.nc`: Ambient temperature data from HERA -- `resources//regions_onshore_base_s_{clusters}.geojson`: Onshore regions -- `resources//dh_areas_base_s_{clusters}.geojson`: District heating areas +- ``data/lake_data/HydroLAKES_polys_v10.gdb/``: Lake polygons from HydroLAKES +- ``data/hera_{year}/ambient_temp_{year}.nc``: Ambient temperature data from HERA +- ``resources//regions_onshore_base_s_{clusters}.geojson``: Onshore regions +- ``resources//dh_areas_base_s_{clusters}.geojson``: District heating areas Outputs ------- -- `resources//heat_source_power_river_water_base_s_{clusters}.csv`: River heating power potentials by region -- `resources//temp_river_water_base_s_{clusters}.nc`: River water temperature profiles by region -- `resources//temp_river_water_base_s_{clusters}_temporal_aggregate.nc`: Temporal aggregated temperature data -- `resources//heat_source_energy_river_water_base_s_{clusters}_temporal_aggregate.nc`: Temporal aggregated energy data +- ``resources//heat_source_power_lake_water_base_s_{clusters}.csv``: Lake heating power potentials by region +- ``resources//temp_lake_water_base_s_{clusters}.nc``: Lake water temperature profiles by region +- ``resources//temp_lake_water_base_s_{clusters}_temporal_aggregate.nc``: Temporal aggregated temperature data +- ``resources//heat_source_energy_lake_water_base_s_{clusters}_temporal_aggregate.nc``: Temporal aggregated energy data """ import gc import logging +import warnings import dask import geopandas as gpd @@ -58,6 +62,12 @@ ) from approximators.lake_water_heat_approximator import LakeWaterHeatApproximator +warnings.filterwarnings( + "ignore", + message=".*organizePolygons\\(\\) received a polygon with more than 100 parts.*", +) + + logger = logging.getLogger(__name__) MEMORY_SAFETY_FACTOR = 0.7 # Use 70% of available memory for Dask arrays @@ -70,24 +80,24 @@ def load_hera_data( miny: float, maxx: float, maxy: float, -) -> dict: +) -> xr.DataArray: """ - Load and concatenate HERA data files for multiple years with spatial clipping. + Load and concatenate HERA ambient temperature data with spatial clipping. Parameters ---------- hera_inputs : dict - Dictionary with year-specific HERA file paths from input_hera_data(). - Expected keys: hera_river_discharge_{year}, hera_ambient_temperature_{year} + Dictionary with year-specific HERA file paths. + Expected keys: hera_ambient_temperature_{year}. snapshots : pd.DatetimeIndex - Target snapshots to select from the combined data + Target snapshots to select from the combined data. minx, miny, maxx, maxy : float - Bounding box coordinates for spatial clipping + Bounding box coordinates for spatial clipping (EPSG:4326). Returns ------- - dict - Dictionary with processed xarray DataArrays for river_discharge and ambient_temperature + xr.DataArray + Ambient temperature data reprojected to EPSG:3035. """ temp_files = [ v for k, v in hera_inputs.items() if k.startswith("hera_ambient_temperature_") @@ -124,54 +134,49 @@ def _create_empty_datasets( snapshots: pd.DatetimeIndex, center_lon: float, center_lat: float ) -> tuple[xr.Dataset, xr.Dataset]: """ - Create empty spatial and temporal aggregate datasets for regions without DH areas. + Create empty datasets for regions without DH areas. When a region has no intersection with district heating areas, we still need - to provide valid datasets with zero values to maintain consistent data structure - across all regions. This prevents errors in downstream processing. + to provide valid datasets with zero values to maintain consistent data structure. Parameters ---------- snapshots : pd.DatetimeIndex - Time snapshots for the spatial aggregate + Time snapshots for the spatial aggregate. center_lon : float - Longitude of region center (for fallback coordinate) + Longitude of region center (for fallback coordinate). center_lat : float - Latitude of region center (for fallback coordinate) + Latitude of region center (for fallback coordinate). Returns ------- tuple[xr.Dataset, xr.Dataset] - Tuple of (spatial_aggregate, temporal_aggregate) datasets with zero values + Tuple of (spatial_aggregate, temporal_aggregate) datasets with zero values. """ - # Create spatial aggregate with time-series of zeros - # This represents no available river heat power/temperature over time spatial_aggregate = xr.Dataset( data_vars={ "total_power": xr.DataArray( - np.zeros(len(snapshots)), # Zero power for all timestamps + np.zeros(len(snapshots)), dims=["time"], coords={"time": snapshots}, ), "average_temperature": xr.DataArray( - np.zeros(len(snapshots)), # Zero temperature for all timestamps + np.zeros(len(snapshots)), dims=["time"], coords={"time": snapshots}, ), } ) - # Create temporal aggregate with single spatial point of zeros - # This represents no available river heat energy/temperature spatially temporal_aggregate = xr.Dataset( data_vars={ "total_energy": xr.DataArray( - [[0.0]], # Zero energy at region center + [[0.0]], dims=["longitude", "latitude"], coords={"longitude": [center_lon], "latitude": [center_lat]}, ), "average_temperature": xr.DataArray( - [[0.0]], # Zero temperature at region center + [[0.0]], dims=["longitude", "latitude"], coords={"longitude": [center_lon], "latitude": [center_lat]}, ), @@ -190,32 +195,35 @@ def get_regional_result( enable_heat_source_maps: bool = False, ) -> dict[str, xr.Dataset]: """ - Calculate river water heat potential for a given region and district heating areas. + Calculate lake water heat potential for a given region. Parameters ---------- hera_inputs : dict Dictionary containing HERA input file paths with year-specific keys. - region : geopandas.GeoSeries + region : gpd.GeoSeries Geographical region for which to compute the heat potential. - dh_areas : geopandas.GeoDataFrame + dh_areas : gpd.GeoDataFrame District heating areas to intersect with the region. + lake_shapes : gpd.GeoDataFrame + Lake polygons from HydroLAKES with Vol_total attribute. snapshots : pd.DatetimeIndex - Time snapshots, used for loading data and for regions without dh_areas - enable_heat_source_maps : bool, default False - Whether to enable heat source mapping. + Time snapshots for temporal alignment. + enable_heat_source_maps : bool, optional + Whether to compute temporal aggregate for heat source maps, + by default False. Returns ------- - dict - Dictionary with keys 'spatial aggregate' and 'temporal aggregate'. - 'spatial aggregate' contains total power and average temperature. - 'temporal aggregate' contains time series for energy and temperature (for analysis/plotting). + dict[str, xr.Dataset] + Dictionary with keys: + - 'spatial aggregate': Dataset with total_power [MW] and + average_temperature [°C] time series. + - 'temporal aggregate' (optional): Dataset with spatial distribution + of total_energy [MWh] and average_temperature [°C]. """ - # Store original region for fallback centroid calculation original_region = region.copy() - # Intersect region with district heating areas intersected_geometry = gpd.overlay( region.to_frame(), dh_areas, @@ -224,95 +232,78 @@ def get_regional_result( region.geometry = intersected_geometry - # Handle empty geometry case (no intersection with DH areas) - # This occurs when a region has no district heating area - # Note: the region could still have district heating, as central heat demand is computed elsewhere + # Return empty datasets for regions without DH area intersection if region.geometry.is_empty.any(): - # Get the center of the original region (before intersection) - # We use the original region to get a meaningful coordinate for the empty datasets - # Project to EPSG:3035 for accurate centroid calculation, then back to EPSG:4326 region_center = ( original_region.to_crs("EPSG:3035").centroid.to_crs("EPSG:4326").iloc[0] ) - center_lon = region_center.x - center_lat = region_center.y - - # Return zero-filled datasets with proper structure spatial_aggregate, temporal_aggregate = _create_empty_datasets( - snapshots, center_lon, center_lat + snapshots, region_center.x, region_center.y ) - return { "spatial aggregate": spatial_aggregate, "temporal aggregate": temporal_aggregate, } - # Process region with valid DH area intersection - # Get bounding box for efficient data clipping minx, miny, maxx, maxy = region.total_bounds - - # Load and concatenate HERA data for all required years with preprocessing ambient_temperature = load_hera_data(hera_inputs, snapshots, minx, miny, maxx, maxy) - # Reproject region to match data CRS for spatial calculations region = region.to_crs("EPSG:3035") lake_shapes = lake_shapes.to_crs("EPSG:3035") lake_water_heat_approximator = LakeWaterHeatApproximator( - ambient_temperature=ambient_temperature, # Air temperature - region=region, # Geographic region of interest - lake_shapes=lake_shapes, # Lake polygons with Vol_total + ambient_temperature=ambient_temperature, + region=region, + lake_shapes=lake_shapes, ) - # Calculate spatial aggregate (time series data for the entire region) - # Contains total_power [MW] and average_temperature [°C] over time spatial_aggregate = lake_water_heat_approximator.get_spatial_aggregate() - # Calculate temporal aggregate only if needed (spatial distribution data for plotting, no time dimension) if enable_heat_source_maps: temporal_aggregate = ( - lake_water_heat_approximator.get_temporal_aggregate().rio.reproject( - "EPSG:4326" - ) # Convert back to WGS84 for output consistency - # .rename({"x": "longitude", "y": "latitude"}) # Standardize coordinate names + lake_water_heat_approximator.get_temporal_aggregate() + .rio.reproject("EPSG:4326") + .compute() ) - temporal_aggregate = temporal_aggregate.compute() else: temporal_aggregate = None - # Compute results immediately to free Dask arrays and enable garbage collection spatial_aggregate = spatial_aggregate.compute() - # Explicitly delete approximator and intermediate arrays - del lake_water_heat_approximator - del ambient_temperature + # Free memory for next region + del lake_water_heat_approximator, ambient_temperature gc.collect() - result = { - "spatial aggregate": spatial_aggregate, - } - - # Only include temporal aggregate if computed + result = {"spatial aggregate": spatial_aggregate} if temporal_aggregate is not None: result["temporal aggregate"] = temporal_aggregate - if enable_heat_source_maps: - del temporal_aggregate return result def set_dask_chunk_size( - n_threads: int, # Number of threads per worker, - memory_mb: int, # Memory per worker in MB - memory_safety_factor=MEMORY_SAFETY_FACTOR, - n_datasets: int = 2, # ambient temperature and river discharge datasets - operation_multiplier: int = 3, # Multiplier for operation overhead + n_threads: int, + memory_mb: int, + memory_safety_factor: float = MEMORY_SAFETY_FACTOR, + n_datasets: int = 1, + operation_multiplier: int = 3, ) -> None: """ - Set the Dask chunk size based on available memory and number of threads. - This function calculates the chunk size for Dask arrays to optimize memory usage - """ + Configure Dask chunk size based on available memory. + Parameters + ---------- + n_threads : int + Number of threads per worker. + memory_mb : int + Memory per worker in MB. + memory_safety_factor : float, optional + Fraction of memory to use, by default 0.7. + n_datasets : int, optional + Number of concurrent datasets in memory, by default 1. + operation_multiplier : int, optional + Multiplier for operation overhead, by default 3. + """ chunk_size = ( memory_mb * memory_safety_factor / n_threads / n_datasets / operation_multiplier ) @@ -324,7 +315,7 @@ def set_dask_chunk_size( from _helpers import mock_snakemake snakemake = mock_snakemake( - "build_river_water_heat_potential", + "build_lake_water_heat_potential", clusters="39", opts="", ll="vopt", From cd58ff8c6da407c36cb3a4c98f2d9663f28994cd Mon Sep 17 00:00:00 2001 From: Amos Schledorn Date: Tue, 23 Dec 2025 17:51:23 +0100 Subject: [PATCH 06/22] feat: include lake water in supply limited heat sources --- scripts/definitions/heat_source.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/scripts/definitions/heat_source.py b/scripts/definitions/heat_source.py index 7ac743c85f..1ce8a9d1e5 100644 --- a/scripts/definitions/heat_source.py +++ b/scripts/definitions/heat_source.py @@ -112,7 +112,11 @@ def source_type(self) -> HeatSourceType: """ if self in [HeatSource.AIR, HeatSource.GROUND, HeatSource.SEA_WATER]: return HeatSourceType.INEXHAUSTIBLE - elif self in [HeatSource.GEOTHERMAL, HeatSource.RIVER_WATER]: + elif self in [ + HeatSource.GEOTHERMAL, + HeatSource.RIVER_WATER, + HeatSource.LAKE_WATER, + ]: return HeatSourceType.SUPPLY_LIMITED elif self == HeatSource.PTES: return HeatSourceType.STORAGE From 39323681c6b7ac75200f48c73cf27445e9d6fdc2 Mon Sep 17 00:00:00 2001 From: Amos Schledorn Date: Tue, 23 Dec 2025 17:51:33 +0100 Subject: [PATCH 07/22] feat: add lake water heat source colors to plotting configuration --- config/plotting.default.yaml | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/config/plotting.default.yaml b/config/plotting.default.yaml index bb31bf9f46..815d051f21 100644 --- a/config/plotting.default.yaml +++ b/config/plotting.default.yaml @@ -470,6 +470,10 @@ plotting: river_water heat preheater: '#4bb9f2' river_water heat utilisation: '#4bb9f2' river_water heat pump: '#4bb9f2' + lake_water heat: '#8db2f7' + lake_water heat preheater: '#8db2f7' + lake_water heat utilisation: '#8db2f7' + lake_water heat pump: '#b1c8f2' sea_water heat: '#0b222e' sea_water heat pump: '#0b222e' ground heat pump: '#2fb537' From 28085a453c985ef624c53e4dcdd45ab940367c2d Mon Sep 17 00:00:00 2001 From: Amos Schledorn Date: Tue, 23 Dec 2025 17:53:35 +0100 Subject: [PATCH 08/22] doc: update docs and release notes --- doc/configtables/sector.csv | 2 +- doc/release_notes.rst | 2 ++ 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/doc/configtables/sector.csv b/doc/configtables/sector.csv index 7a945b9702..0e8c932d72 100644 --- a/doc/configtables/sector.csv +++ b/doc/configtables/sector.csv @@ -59,7 +59,7 @@ district_heating,--,, -- -- constant_temperature_celsius,C,float,Constant temperature assumed for geothermal heat sources. Used for scaling heat potentials based on actual vs assumed temperature differences. -- fallback_ptx_heat_losses,--,float,"Non-recoverable heat loss factor for PTX processes without efficiency-heat data in technology-data. heat_sources,--,, --- urban central,--,"List of heat sources for heat pumps in urban central heating. Must be one of [air, geothermal, ptes, river_water, sea_water, electrolysis_excess, fuel_cell_excess, fischer_tropsch_excess, haber_bosch_excess, sabatier_excess, methanolisation_excess]", +-- urban central,--,"List of heat sources for heat pumps in urban central heating. Must be one of [air, geothermal, ptes, river_water, lake_water, sea_water, electrolysis_excess, fuel_cell_excess, fischer_tropsch_excess, haber_bosch_excess, sabatier_excess, methanolisation_excess]", -- urban decentral,--,"List of heat sources for heat pumps in urban decentral heating. Must be one of [air]", -- rural,--,"List of heat sources for heat pumps in rural heating. Must be one of [air, ground]", heat_source_temperatures,--,,Dictionary of temperatures (°C) for each heat source. Used for computing source utilisation. diff --git a/doc/release_notes.rst b/doc/release_notes.rst index 5a71551fac..3b6015c1b1 100644 --- a/doc/release_notes.rst +++ b/doc/release_notes.rst @@ -9,6 +9,8 @@ Release Notes Upcoming Release ================ +* Feat: added support for lake-water sourced heat pumps in district heating, based on data from HydroLakes and a methods from Triebs (https://github.com/PyPSA/pypsa-eur/pull/1951) + * Refactor: Integrated excess heat from Power-to-X processes into the new heat-source structure and moved some code from `scripts/prepare_sector_network.py` to `scripts/def/heat_source.py`. Also updated PtX excess heat efficiencies (https://github.com/PyPSA/pypsa-eur/pull/1944). * Update heat source handling in `prepare_sector_network` and introduce preheating of heat sources for more realistic system integrations (https://github.com/PyPSA/pypsa-eur/pull/1893). From b48210fb68c8bc50a0f42da7e9255a11cdfbaae8 Mon Sep 17 00:00:00 2001 From: Amos Schledorn Date: Tue, 23 Dec 2025 17:54:29 +0100 Subject: [PATCH 09/22] feat: add lake-water HPs to test configs --- config/test/config.myopic.yaml | 1 + config/test/config.overnight.yaml | 1 + config/test/config.perfect.yaml | 1 + 3 files changed, 3 insertions(+) diff --git a/config/test/config.myopic.yaml b/config/test/config.myopic.yaml index 4f8e639e74..0bb1d746d9 100644 --- a/config/test/config.myopic.yaml +++ b/config/test/config.myopic.yaml @@ -40,6 +40,7 @@ sector: - air - geothermal - river_water + - lake_water - sea_water - ptes urban decentral: diff --git a/config/test/config.overnight.yaml b/config/test/config.overnight.yaml index 936cc2afe8..3eef0cd8f6 100644 --- a/config/test/config.overnight.yaml +++ b/config/test/config.overnight.yaml @@ -80,6 +80,7 @@ sector: - air - geothermal - river_water + - lake_water - sea_water - ptes urban decentral: diff --git a/config/test/config.perfect.yaml b/config/test/config.perfect.yaml index 92cbc2535e..cae7bc8790 100644 --- a/config/test/config.perfect.yaml +++ b/config/test/config.perfect.yaml @@ -52,6 +52,7 @@ sector: - ptes - river_water - sea_water + - lake_water - geothermal - ptes urban decentral: From 767371b3dbcb15b1571c26626389f2af2237e57a Mon Sep 17 00:00:00 2001 From: Amos Schledorn Date: Tue, 6 Jan 2026 10:38:17 +0100 Subject: [PATCH 10/22] integrate HydroLAKES data into new data structure --- data/versions.csv | 2 +- doc/configtables/data.csv | 5 ++++- 2 files changed, 5 insertions(+), 2 deletions(-) diff --git a/data/versions.csv b/data/versions.csv index 4d4f7ffc2a..80625356fc 100644 --- a/data/versions.csv +++ b/data/versions.csv @@ -100,4 +100,4 @@ "jrc_ardeco","primary","2021","['latest', 'supported']","https://territorial.ec.europa.eu/ardeco-api-v2/rest/export/","" "jrc_ardeco","archive","2021","['latest', 'supported']","https://zenodo.org/records/17249457/files","" "lake_data","primary","v10","['latest', 'supported']","https://data.hydrosheds.org/file/hydrolakes/HydroLAKES_polys_v10.gdb.zip","" -"lake_data","archive","v10","['latest', 'supported']","https://data.hydrosheds.org/file/hydrolakes/HydroLAKES_polys_v10.gdb.zip","" \ No newline at end of file +"lake_data","archive","v10","['latest', 'supported']","https://zenodo.org/records/18153550/files/HydroLAKES_polys_v10.gdb.zip?download=1","" \ No newline at end of file diff --git a/doc/configtables/data.csv b/doc/configtables/data.csv index 2f0b075671..40cad6b490 100644 --- a/doc/configtables/data.csv +++ b/doc/configtables/data.csv @@ -136,4 +136,7 @@ geothermal_heat_utilisation_potentials,,, -- version,str,,"Version of the Geothermal Heat Utilisation Potentials" jrc_ardeco,,, -- source,str,"{archive, primary}","Source of annual regional database for EU regions" --- version,str,,"Version of the JRC ARDECO dataset" \ No newline at end of file +-- version,str,,"Version of the JRC ARDECO dataset" +lake_data,,, +-- source,str,"{archive, primary}","Source of global lake polygons" +-- version,str,,"Version of the HydroLAKES dataset" \ No newline at end of file From aa8bcee53aa159fc6337b6d30c81791c447a16fc Mon Sep 17 00:00:00 2001 From: Amos Schledorn <60692940+amos-schledorn@users.noreply.github.com> Date: Tue, 6 Jan 2026 10:45:07 +0100 Subject: [PATCH 11/22] Update config/config.default.yaml Co-authored-by: Fabian Neumann --- config/config.default.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/config/config.default.yaml b/config/config.default.yaml index 398445572e..27d107e249 100644 --- a/config/config.default.yaml +++ b/config/config.default.yaml @@ -1150,7 +1150,7 @@ data: source: archive version: latest lake_data: - source: archive + source: primary version: latest secrets: From 18656fbf5e696421f91b151358ff9169ebb5c22d Mon Sep 17 00:00:00 2001 From: Amos Schledorn Date: Tue, 6 Jan 2026 14:30:00 +0100 Subject: [PATCH 12/22] feat: add rule to retrieve lake data from HydroLAKES for Belgium --- rules/retrieve.smk | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/rules/retrieve.smk b/rules/retrieve.smk index 09d4936332..0ea0d4ba25 100755 --- a/rules/retrieve.smk +++ b/rules/retrieve.smk @@ -1179,3 +1179,20 @@ if (LAKE_DATA_DATASET := dataset_version("lake_data"))["source"] in [ output["zip_file"], LAKE_DATA_DATASET["folder"], ) + + rule retrieve_lake_data_test_cutout: + input: + zip_file=storage( + f"https://zenodo.org/records/18163260/files/HydroLAKES_polys_v10_belgium.gdb.zip?download=1" + ), + output: + zip_file=f"{LAKE_DATA_DATASET['folder']}/HydroLAKES_polys_v10_belgium.gdb.zip", + lake_data=directory( + f"{LAKE_DATA_DATASET['folder']}/HydroLAKES_polys_v10_belgium.gdb" + ), + run: + copy2(input["zip_file"], output["zip_file"]) + unpack_archive( + output["zip_file"], + LAKE_DATA_DATASET["folder"], + ) From 0a89f7c139a0ace1bfd3955620cf522a55f7745e Mon Sep 17 00:00:00 2001 From: Amos Schledorn Date: Wed, 7 Jan 2026 10:53:18 +0100 Subject: [PATCH 13/22] fall back to Europe lake data as archive --- config/config.default.yaml | 2 +- data/versions.csv | 2 +- rules/retrieve.smk | 17 ----------------- 3 files changed, 2 insertions(+), 19 deletions(-) diff --git a/config/config.default.yaml b/config/config.default.yaml index 27d107e249..398445572e 100644 --- a/config/config.default.yaml +++ b/config/config.default.yaml @@ -1150,7 +1150,7 @@ data: source: archive version: latest lake_data: - source: primary + source: archive version: latest secrets: diff --git a/data/versions.csv b/data/versions.csv index 80625356fc..c4ea64ac89 100644 --- a/data/versions.csv +++ b/data/versions.csv @@ -100,4 +100,4 @@ "jrc_ardeco","primary","2021","['latest', 'supported']","https://territorial.ec.europa.eu/ardeco-api-v2/rest/export/","" "jrc_ardeco","archive","2021","['latest', 'supported']","https://zenodo.org/records/17249457/files","" "lake_data","primary","v10","['latest', 'supported']","https://data.hydrosheds.org/file/hydrolakes/HydroLAKES_polys_v10.gdb.zip","" -"lake_data","archive","v10","['latest', 'supported']","https://zenodo.org/records/18153550/files/HydroLAKES_polys_v10.gdb.zip?download=1","" \ No newline at end of file +"lake_data","archive","v10","['latest', 'supported']","https://zenodo.org/records/18164492/files/HydroLAKES_polys_v10_europe.gdb.zip?download=1","" \ No newline at end of file diff --git a/rules/retrieve.smk b/rules/retrieve.smk index 0ea0d4ba25..09d4936332 100755 --- a/rules/retrieve.smk +++ b/rules/retrieve.smk @@ -1179,20 +1179,3 @@ if (LAKE_DATA_DATASET := dataset_version("lake_data"))["source"] in [ output["zip_file"], LAKE_DATA_DATASET["folder"], ) - - rule retrieve_lake_data_test_cutout: - input: - zip_file=storage( - f"https://zenodo.org/records/18163260/files/HydroLAKES_polys_v10_belgium.gdb.zip?download=1" - ), - output: - zip_file=f"{LAKE_DATA_DATASET['folder']}/HydroLAKES_polys_v10_belgium.gdb.zip", - lake_data=directory( - f"{LAKE_DATA_DATASET['folder']}/HydroLAKES_polys_v10_belgium.gdb" - ), - run: - copy2(input["zip_file"], output["zip_file"]) - unpack_archive( - output["zip_file"], - LAKE_DATA_DATASET["folder"], - ) From e3132d6ebc5a933983b34ebefffe594ea779626d Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Wed, 7 Jan 2026 10:58:43 +0000 Subject: [PATCH 14/22] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- .../approximators/lake_water_heat_approximator.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/scripts/build_surface_water_heat_potentials/approximators/lake_water_heat_approximator.py b/scripts/build_surface_water_heat_potentials/approximators/lake_water_heat_approximator.py index 9ce85eae7a..25840531c1 100644 --- a/scripts/build_surface_water_heat_potentials/approximators/lake_water_heat_approximator.py +++ b/scripts/build_surface_water_heat_potentials/approximators/lake_water_heat_approximator.py @@ -4,7 +4,6 @@ """Lake water heat approximator for district heating systems.""" from functools import cached_property -from typing import Union import geopandas as gpd import shapely @@ -53,7 +52,7 @@ class LakeWaterHeatApproximator(SurfaceWaterHeatApproximator): def __init__( self, ambient_temperature: xr.DataArray, - region: Union[shapely.geometry.polygon.Polygon, gpd.GeoSeries], + region: shapely.geometry.polygon.Polygon | gpd.GeoSeries, lake_shapes: gpd.GeoDataFrame, delta_t_max: float = 1, min_outlet_temperature: float = 1, From f7003763441ce33ee6e9137d0521957da6918271 Mon Sep 17 00:00:00 2001 From: Amos Schledorn Date: Wed, 7 Jan 2026 12:00:38 +0100 Subject: [PATCH 15/22] fix: clarify preheating logic in get_sink_inlet_temperature function --- scripts/build_cop_profiles/run.py | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/scripts/build_cop_profiles/run.py b/scripts/build_cop_profiles/run.py index 076478e320..949033c564 100644 --- a/scripts/build_cop_profiles/run.py +++ b/scripts/build_cop_profiles/run.py @@ -237,6 +237,9 @@ def get_sink_inlet_temperature( temperature. When preheating is not used (source <= return), the heat pump receives water at return temperature and heats it to forward temperature. + When source temperature > return temperature, preheater is used: preheater raises return flow, heat pump inlet is at forward temp. + When source temperature <= return temperature, no preheating: heat pump inlet is at return temperature. + Parameters ---------- heat_source_name : str @@ -255,10 +258,6 @@ def get_sink_inlet_temperature( """ heat_source = HeatSource(heat_source_name) if heat_source.requires_preheater: - # When source temperature > return temperature, preheater is used: - # preheater raises return flow, heat pump inlet is at forward temp. - # When source temperature <= return temperature, no preheating: - # heat pump inlet is at return temperature. return central_heating_forward_temperature.where( central_heating_return_temperature < source_temperature, central_heating_return_temperature, @@ -320,7 +319,7 @@ def get_sink_inlet_temperature( heat_source=heat_source_name, source_inlet_temperature_celsius=source_inlet_temperature_celsius, sink_outlet_temperature_celsius=central_heating_forward_temperature, - sink_inlet_temperature_celsius=central_heating_return_temperature, + sink_inlet_temperature_celsius=sink_inlet_temperature_celsius, ) cop_this_system_type.append(cop_da) cop_all_system_types.append( From 1ccf297eafa018b52f99fd60c89a7f4514c9f1c9 Mon Sep 17 00:00:00 2001 From: Amos Schledorn Date: Tue, 3 Feb 2026 17:50:18 +0100 Subject: [PATCH 16/22] feat: add lake data source configuration and heat source temperature settings --- data/versions.csv | 2 ++ scripts/lib/validation/config/data.py | 4 ++++ scripts/lib/validation/config/sector.py | 16 ++++++++++++++++ 3 files changed, 22 insertions(+) diff --git a/data/versions.csv b/data/versions.csv index 94eabacb14..cab97fae35 100644 --- a/data/versions.csv +++ b/data/versions.csv @@ -137,3 +137,5 @@ wdpa_marine,unknown,primary,latest supported,2025-12-02,"WDPA maritime changes i wdpa_marine,2025-07,archive,latest supported,2025-12-02,"We are legally not allowed to redistribute this dataset, luckily the web archive is keeping copies of it.",https://web.archive.org/web/20250715084308if_/https%3A//d1gam3xoknrgr2.cloudfront.net/current/WDPA_WDOECM_Jul2025_Public_marine_shp.zip worldbank_urban_population,unknown,primary,latest might-work,2025-12-02,"This is the original World Bank API link, which is sometimes updated; it is not guaranteed to work with the current codebase and data changes without notice.",https://api.worldbank.org/v2/en/indicator/SP.URB.TOTL.IN.ZS?downloadformat=csv worldbank_urban_population,2025-08-14,archive,latest supported,2026-01-13,,https://data.pypsa.org/workflows/eur/worldbank_urban_population/2025-08-14/API_SP.URB.TOTL.IN.ZS_DS2_en_csv_v2_22447.zip +lake_data,v10,primary,latest supported,2025-12-02,,https://data.hydrosheds.org/file/hydrolakes/HydroLAKES_polys_v10.gdb.zip +lake_data,v10,archive,latest supported,2026-01-13,,https://zenodo.org/records/18164492/files/HydroLAKES_polys_v10_europe.gdb.zip?download=1 diff --git a/scripts/lib/validation/config/data.py b/scripts/lib/validation/config/data.py index 3b814f8500..b7d03f59ac 100644 --- a/scripts/lib/validation/config/data.py +++ b/scripts/lib/validation/config/data.py @@ -243,3 +243,7 @@ class DataConfig(BaseModel): default_factory=_DataSourceConfig, description="Entsoepy bidding zones data source configuration.", ) + lake_data: _DataSourceConfig = Field( + default_factory=_DataSourceConfig, + description="Lake data source configuration.", + ) diff --git a/scripts/lib/validation/config/sector.py b/scripts/lib/validation/config/sector.py index b69dfd2b2f..74bc678224 100644 --- a/scripts/lib/validation/config/sector.py +++ b/scripts/lib/validation/config/sector.py @@ -185,6 +185,22 @@ class _DistrictHeatingConfig(ConfigModel): default_factory=lambda: {"buffer": 1000, "handle_missing_countries": "fill"}, description="District heating areas settings.", ) + heat_source_temperatures: dict[str, float] = Field( + default_factory=lambda: { + "geothermal": 65, + "electrolysis_waste": 70, + "fuel_cell_waste": 70, + "fischer_tropsch_waste": 200, + "haber_bosch_waste": 200, + "sabatier_waste": 200, + "methanolisation_waste": 200, + }, + description="Temperature in °C for each heat source. Determines whether heat can be used directly (T_source > T_forward), for preheating (T_return < T_source < T_forward), or requires boosting via heat pumps.", + ) + fallback_ptx_heat_losses: float = Field( + 0.05, + description="Default heat loss fraction for PtX processes when waste heat recovery is enabled but no specific loss value is provided.", + ) geothermal: dict[str, Any] = Field( default_factory=lambda: {"constant_temperature_celsius": 65}, description=( From 16450e49fd96290d2d4ade2dceafe0df2ede19af Mon Sep 17 00:00:00 2001 From: Amos Schledorn Date: Mon, 23 Mar 2026 15:25:22 +0100 Subject: [PATCH 17/22] fix lake data retrieval destination folder --- rules/retrieve.smk | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/rules/retrieve.smk b/rules/retrieve.smk index 4ac860cf4b..bd9e391482 100755 --- a/rules/retrieve.smk +++ b/rules/retrieve.smk @@ -1708,5 +1708,5 @@ if (LAKE_DATA_DATASET := dataset_version("lake_data"))["source"] in [ copy2(input["zip_file"], output["zip_file"]) unpack_archive( output["zip_file"], - LAKE_DATA_DATASET["folder"], + output["lake_data"], ) From f7af50ff03f21e5745eac627f95faedcfe1fb13f Mon Sep 17 00:00:00 2001 From: Amos Schledorn Date: Tue, 24 Mar 2026 12:58:23 +0100 Subject: [PATCH 18/22] update lake data retrieval --- config/config.default.yaml | 2 +- rules/retrieve.smk | 1 + 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/config/config.default.yaml b/config/config.default.yaml index 936330feec..a0fef79a73 100644 --- a/config/config.default.yaml +++ b/config/config.default.yaml @@ -1456,7 +1456,7 @@ data: source: archive version: latest lake_data: - source: archive + source: primary version: latest # docs in https://pypsa-eur.readthedocs.io/en/latest/configuration.html#overpass_api diff --git a/rules/retrieve.smk b/rules/retrieve.smk index bd9e391482..44f2a247e0 100755 --- a/rules/retrieve.smk +++ b/rules/retrieve.smk @@ -1703,6 +1703,7 @@ if (LAKE_DATA_DATASET := dataset_version("lake_data"))["source"] in [ zip_file=f"{LAKE_DATA_DATASET['folder']}/HydroLAKES_polys_v10.gdb.zip", lake_data=directory( f"{LAKE_DATA_DATASET['folder']}/HydroLAKES_polys_v10.gdb" + f"{LAKE_DATA_DATASET['folder']}/HydroLAKES_polys_v10.gdb/HydroLAKES_polys_v10.gdb" ), run: copy2(input["zip_file"], output["zip_file"]) From cec3d7820d95f73d51c95d561adb1fe4eb0b280c Mon Sep 17 00:00:00 2001 From: Amos Schledorn Date: Tue, 24 Mar 2026 12:59:19 +0100 Subject: [PATCH 19/22] fix: correct lake data output path in retrieval script --- rules/retrieve.smk | 1 - 1 file changed, 1 deletion(-) diff --git a/rules/retrieve.smk b/rules/retrieve.smk index 44f2a247e0..2c1976854f 100755 --- a/rules/retrieve.smk +++ b/rules/retrieve.smk @@ -1702,7 +1702,6 @@ if (LAKE_DATA_DATASET := dataset_version("lake_data"))["source"] in [ output: zip_file=f"{LAKE_DATA_DATASET['folder']}/HydroLAKES_polys_v10.gdb.zip", lake_data=directory( - f"{LAKE_DATA_DATASET['folder']}/HydroLAKES_polys_v10.gdb" f"{LAKE_DATA_DATASET['folder']}/HydroLAKES_polys_v10.gdb/HydroLAKES_polys_v10.gdb" ), run: From 638cd58c49897043725455906677b92d04cc6480 Mon Sep 17 00:00:00 2001 From: Amos Schledorn Date: Wed, 25 Mar 2026 15:02:21 +0100 Subject: [PATCH 20/22] fix handling of small regions for lake heat potentials --- .../lake_water_heat_approximator.py | 25 ++++++++++++++++--- .../build_lake_water_heat_potential.py | 4 +-- 2 files changed, 24 insertions(+), 5 deletions(-) diff --git a/scripts/build_surface_water_heat_potentials/approximators/lake_water_heat_approximator.py b/scripts/build_surface_water_heat_potentials/approximators/lake_water_heat_approximator.py index 25840531c1..de15bd2b55 100644 --- a/scripts/build_surface_water_heat_potentials/approximators/lake_water_heat_approximator.py +++ b/scripts/build_surface_water_heat_potentials/approximators/lake_water_heat_approximator.py @@ -8,6 +8,9 @@ import geopandas as gpd import shapely import xarray as xr +import numpy as np +import rioxarray +import logging from scripts.build_surface_water_heat_potentials.approximators.river_water_heat_approximator import ( RiverWaterHeatApproximator, @@ -16,6 +19,7 @@ SurfaceWaterHeatApproximator, ) +logger = logging.getLogger(__name__) class LakeWaterHeatApproximator(SurfaceWaterHeatApproximator): """ @@ -184,6 +188,9 @@ def _air_temperature_in_lakes( """ Clip ambient temperature raster to lake areas. + Falls back to nearest-neighbor sampling at the centroid if the + geometries are too small to contain any raster cell center. + Parameters ---------- eligible_lake_parts : gpd.GeoDataFrame @@ -194,9 +201,21 @@ def _air_temperature_in_lakes( xr.DataArray Ambient temperature clipped to lake areas. """ - return self.ambient_temperature.rio.clip( - eligible_lake_parts.geometry, eligible_lake_parts.crs, drop=True - ) + if eligible_lake_parts.empty: + logger.warning( + "No raster cells found within lake geometries. Returning NaN." + ) + return xr.full_like(self.ambient_temperature.sel( + x=[self.region.centroid.x[0]], y=[self.region.centroid.y[0]], method="nearest" + ), fill_value=np.nan) + else: + return self.ambient_temperature.rio.clip( + eligible_lake_parts.geometry.buffer(self._data_resolution), # Buffer to ensure we capture nearby raster cells in case there is no raster cell center within the lake polygon + eligible_lake_parts.crs, + drop=True + ) + + @cached_property def _water_temperature_in_region_raster(self) -> xr.DataArray: diff --git a/scripts/build_surface_water_heat_potentials/build_lake_water_heat_potential.py b/scripts/build_surface_water_heat_potentials/build_lake_water_heat_potential.py index bc2ca1a8b2..c6f1f3aadb 100644 --- a/scripts/build_surface_water_heat_potentials/build_lake_water_heat_potential.py +++ b/scripts/build_surface_water_heat_potentials/build_lake_water_heat_potential.py @@ -356,6 +356,7 @@ def set_dask_chunk_size( # Process regions sequentially but with multi-threaded Dask operations results = [] for i, region_name in enumerate(regions_onshore.index, 1): + logger.info(f"Processing region {i}/{len(regions_onshore)}: {region_name}") # Extract region geometry and create a copy to avoid modification conflicts region = gpd.GeoSeries(regions_onshore.loc[region_name].copy(deep=True)) @@ -381,7 +382,7 @@ def set_dask_chunk_size( region_name: res["spatial aggregate"]["total_power"].to_pandas() for region_name, res in zip(regions_onshore.index, results) } - ).dropna() + ).fillna(0) power = power.reindex( snapshots, method="nearest" @@ -404,7 +405,6 @@ def set_dask_chunk_size( dim="name", ) .assign_coords(name=regions_onshore.index) - .dropna(dim="time") ) # Align temperature data to snapshots, use nearest to handle any minor decimal differences From db9e37c92f54679ffa14dabc832e57bcf26b9b77 Mon Sep 17 00:00:00 2001 From: Amos Schledorn <60692940+amos-schledorn@users.noreply.github.com> Date: Fri, 12 Jun 2026 14:19:39 +0200 Subject: [PATCH 21/22] fix flag for ignoring missing geothermal data in non-EU countries (#2127) * fix flag for ignoring missing geothermal data in non-EU countries * update lake data retrieval --- config/config.default.yaml | 1 + rules/build_sector.smk | 6 ++---- scripts/build_geothermal_heat_potential.py | 8 ++++---- scripts/lib/validation/config/sector.py | 5 +++++ 4 files changed, 12 insertions(+), 8 deletions(-) diff --git a/config/config.default.yaml b/config/config.default.yaml index a0fef79a73..b49726f95a 100644 --- a/config/config.default.yaml +++ b/config/config.default.yaml @@ -696,6 +696,7 @@ sector: sabatier_waste: 200 methanolisation_waste: 200 fallback_ptx_heat_losses: 0.05 + ignore_missing_geothermal_data: false heat_sources: urban central: - air diff --git a/rules/build_sector.smk b/rules/build_sector.smk index ccad6b91b0..2c57ce9490 100755 --- a/rules/build_sector.smk +++ b/rules/build_sector.smk @@ -329,12 +329,10 @@ rule build_geothermal_heat_potential: "heat_source_temperatures", "geothermal", ), - ignore_missing_regions=config_provider( + ignore_missing_geothermal_data=config_provider( "sector", "district_heating", - "limited_heat_sources", - "geothermal", - "ignore_missing_regions", + "ignore_missing_geothermal_data", ), heat_source_cooling=config_provider( "sector", "district_heating", "heat_source_cooling" diff --git a/scripts/build_geothermal_heat_potential.py b/scripts/build_geothermal_heat_potential.py index c5d20fbdfa..e7c02d0f49 100644 --- a/scripts/build_geothermal_heat_potential.py +++ b/scripts/build_geothermal_heat_potential.py @@ -242,8 +242,8 @@ def get_heat_source_power( supply_potentials: gpd.GeoDataFrame, full_load_hours: float, input_unit: str, + ignore_missing_data: bool, output_unit: str = "MWh", - ignore_missing_regions: bool = False, ) -> pd.DataFrame: """ Get the heat source power from supply potentials. @@ -309,7 +309,7 @@ def get_heat_source_power( if not non_covered_regions.empty: if all(non_covered_regions.str.contains("|".join(not_eu_27))): - if ignore_missing_regions: + if ignore_missing_data: logger.warning( f"The onshore regions outside EU 27 ({non_covered_regions.to_list()}) have no heat source power. Filling with zeros." ) @@ -318,7 +318,7 @@ def get_heat_source_power( ) else: raise ValueError( - f"The onshore regions outside EU 27 {non_covered_regions.to_list()} have no heat source power. Set the ignore_missing_regions parameter in the config to true if you want to include these countries in your analysis despite missing geothermal data." + f"The onshore regions outside EU 27 {non_covered_regions.to_list()} have no heat source power. Set the ignore_missing_geothermal_data parameter in the config to true if you want to include these countries in your analysis despite missing geothermal data." ) else: raise ValueError( @@ -390,7 +390,7 @@ def get_heat_source_power( supply_potentials=geothermal_supply_potentials, full_load_hours=FULL_LOAD_HOURS, input_unit=input_unit, - ignore_missing_regions=snakemake.params.ignore_missing_regions, + ignore_missing_data=snakemake.params.ignore_missing_geothermal_data, ) # Load temperature profiles for scaling diff --git a/scripts/lib/validation/config/sector.py b/scripts/lib/validation/config/sector.py index 1ca6f49d04..68f337978f 100644 --- a/scripts/lib/validation/config/sector.py +++ b/scripts/lib/validation/config/sector.py @@ -189,6 +189,11 @@ class _DistrictHeatingConfig(ConfigModel): 0.05, description="Default heat loss fraction for PtX processes when waste heat recovery is enabled but no specific loss value is provided.", ) + ignore_missing_geothermal_data: bool = Field( + False, + description="If true, missing geothermal data for non-EU countries will be ignored. If false, an error will be raised if countries with missing data are modelled while geothermal heat is included as a heat source.", + ) + heat_source_temperatures: dict[str, float] = Field( default_factory=lambda: { "geothermal": 65, From cddc856b895eaeb1e6ac8b048e570d38fa2134ff Mon Sep 17 00:00:00 2001 From: Amos Schledorn Date: Wed, 17 Jun 2026 11:32:14 +0200 Subject: [PATCH 22/22] fix: correct casing for heat source waste types in HeatSystem enum --- scripts/definitions/heat_system.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/scripts/definitions/heat_system.py b/scripts/definitions/heat_system.py index 76b50285da..67ffcc392b 100644 --- a/scripts/definitions/heat_system.py +++ b/scripts/definitions/heat_system.py @@ -235,12 +235,12 @@ def heat_pump_costs_name(self, heat_source: HeatSource | str) -> str: HeatSource.SEA_WATER, HeatSource.RIVER_WATER, HeatSource.LAKE_WATER, - HeatSource.ELECTROLYSIS_waste, - HeatSource.FISCHER_TROPSCH_waste, - HeatSource.SABATIER_waste, - HeatSource.HABER_BOSCH_waste, - HeatSource.METHANOLISATION_waste, - HeatSource.FUEL_CELL_waste, + HeatSource.ELECTROLYSIS_WASTE, + HeatSource.FISCHER_TROPSCH_WASTE, + HeatSource.SABATIER_WASTE, + HeatSource.HABER_BOSCH_WASTE, + HeatSource.METHANOLISATION_WASTE, + HeatSource.FUEL_CELL_WASTE, ]: return f"{self.central_or_decentral} excess-heat-sourced heat pump" else: