Use Global Buildings datasets to estimate solar-rooftop potentials (#1629)

* Use global buildings dataset to estimate solar rooftop potentials

* add to test runs in config.sector.yaml

* split workflow by countries to reduce memory

* simplify buildings during download, reduce memory by 80%

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* add a simplified model to test runs, add release notes

* use parquet to speed up savings and retrieval, add documentations in scripts

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* Add a definition for install_ration

* Add an explanation for tolerance

Author: virio-andreyana

Snakefile

@@ -163,6 +163,17 @@ if config["enable"].get("retrieve_databundle", True):
             "scripts/retrieve_databundle_light.py"
 
 
+if config["enable"].get("download_global_buildings", True):
+
+    rule download_global_buildings:
+        params:
+            crs=config["crs"],
+        output:
+            "data/global_buildings/{country}_global_buildings_raw.parquet",
+        script:
+            "scripts/download_global_buildings.py"
+
+
 if config["enable"].get("download_osm_data", True):
 
     rule download_osm_data:
@@ -695,6 +706,41 @@ rule cluster_network:
         "scripts/cluster_network.py"
 
 
+solar_rooftop_config = config["sector"]["solar_rooftop"]
+if isinstance(solar_rooftop_config, dict):
+    solar_rooftop_enable = (
+        solar_rooftop_config["enable"] and solar_rooftop_config["use_building_size"]
+    )
+    solar_rooftop_params = {
+        "solar_rooftop_enable": solar_rooftop_enable,
+        "install_ratio": solar_rooftop_config["install_ratio"],
+        "tolerance": solar_rooftop_config["tolerance"],
+    }
+else:
+    solar_rooftop_params = {}
+    solar_rooftop_enable = config["sector"]["solar_rooftop"]
+
+
+rule cluster_global_buildings:
+    params:
+        **solar_rooftop_params,
+        crs=config["crs"],
+    input:
+        country_buildings="data/global_buildings/{country}_global_buildings_raw.parquet",
+        regions_onshore="resources/"
+        + RDIR
+        + "bus_regions/regions_onshore_elec_s{simpl}_{clusters}.geojson",
+    output:
+        solar_rooftop_layout=branch(
+            solar_rooftop_enable,
+            "resources/"
+            + RDIR
+            + "solar_rooftop/solar_rooftop_layout_elec_s{simpl}_{clusters}_{country}.csv",
+        ),
+    script:
+        "scripts/cluster_global_buildings.py"
+
+
 if config["augmented_line_connection"].get("add_to_snakefile") == True:
 
     rule augmented_line_connections:
@@ -1109,6 +1155,16 @@ rule prepare_sector_network:
     input:
         **branch(sector_enable["land_transport"], TRANSPORT),
         **branch(sector_enable["heat"], HEAT),
+        **branch(
+            solar_rooftop_enable,
+            {
+                f"solar_rooftop_layout_{country}": "resources/"
+                + RDIR
+                + "solar_rooftop/solar_rooftop_layout_elec_s{simpl}_{clusters}_"
+                + f"{country}.csv"
+                for country in config["countries"]
+            },
+        ),
         network=RESDIR
         + "prenetworks/elec_s{simpl}_{clusters}_ec_l{ll}_{opts}_{sopts}_{planning_horizons}_{discountrate}_{demand}_presec.nc",
         costs="resources/" + RDIR + "costs_{planning_horizons}.csv",

config.default.yaml

@@ -23,6 +23,7 @@ enable:
   retrieve_databundle_sector: true
   retrieve_cost_data: true # true: retrieves cost data from technology data and saves in resources/costs.csv, false: uses cost data in data/costs.csv
   download_osm_data: true # If 'true', OpenStreetMap data will be downloaded for the above given countries
+  download_global_buildings: false # If 'true', GlobalMLBuildingFootprints data will be downloaded for the above given countries
   build_natura_raster: false # If 'true', then an exclusion raster will be build. Otherwise use pregenerated raster.
   build_cutout: false
   # If "build_cutout" : true, then environmental data is extracted according to `snapshots` date range and `countries`
@@ -661,7 +662,22 @@ sector:
   efficiency_heat_gas_to_elec: 0.9
 
   electricity_distribution_grid: true # adds low voltage buses and shifts AC loads, BEVs, heat pumps, and resistive heaters, micro CHPs to low voltage buses if technologies are present
-  solar_rooftop: true # adds distribution side customer rooftop PV (only work if electricity_distribution_grid: true)
+  solar_rooftop: # adds distribution side customer rooftop PV (only work if electricity_distribution_grid: true)
+    enable: true
+    kW_per_m2: 0.1
+    m2_per_person: 20
+    use_building_size: false
+    # proportion of rooftop suitable for PV installation
+    install_ratio:
+      0: 0
+      10: 0.3
+      100: 0.36
+      200: 0.41
+      450: 0.49
+      2300: 0.66
+     # maximum distance [km] to allocate a building within the nearest bus shapes
+    tolerance: 100
+
   home_battery: true # adds home batteries to low voltage buses ((only work if electricity_distribution_grid: true)
   transmission_efficiency:
     electricity distribution grid:

configs/bundle_config.yaml

@@ -46,6 +46,18 @@ databundles:
     - data/gadm/gadm36_NGA/gadm36_NGA.gpkg  # needed in sector-coupled model
     - data/gadm/gadm36_BEN/gadm36_BEN.gpkg  # needed in sector-coupled model
 
+  # tutorial bundle specific for Nigeria and Benin only (TODO: merge with above)
+  bundle_tutorial_NGBJ_global_buildings:
+    countries: [NG, BJ]
+    tutorial: true
+    category: global_buildings
+    destination: "data"
+    urls:
+      gdrive: https://drive.google.com/file/d/1nnseF1A740Gp0IaryeCQ8IRcEF91-aUe/view?usp=drive_link
+    output:
+    - data/global_buildings/NG_global_buildings_raw.parquet  # needed in cluster_global_buildings
+    - data/global_buildings/BJ_global_buildings_raw.parquet  # needed in cluster_global_buildings
+
   # tutorial bundle specific for Botswana only
   bundle_tutorial_BW:
     countries: [BW]

doc/release_notes.rst

@@ -17,6 +17,8 @@ This part of documentation collects descriptive release notes to capture the mai
 
 * Added a hot-fix to handle UNSD data downtime causing CI to fail `PR #1653 <https://github.com/pypsa-meets-earth/pypsa-earth/pull/1653>`__
 
+* Add an option to estimate solar-rooftop potentials using `GlobalMLBuildingFootprints <https://github.com/microsoft/GlobalMLBuildingFootprints>`_, a simplified data is provided for tutorials, `PR #1629 <https://github.com/pypsa-meets-earth/pypsa-earth/pull/1629>`__
+
 * Align PyPSA-Earth costs with the reference units used in `technology-data`, avoiding discrepancies when combining technologies with different original currency years `PR #1604 <https://github.com/pypsa-meets-earth/pypsa-earth/pull/1604>`__
 
 * Add an option to redefine countries into subregions in ``cluster_networks`` `PR #1542 <https://github.com/pypsa-meets-earth/pypsa-earth/pull/1542>`__

scripts/cluster_global_buildings.py

@@ -0,0 +1,129 @@
+# -*- coding: utf-8 -*-
+# SPDX-FileCopyrightText:  PyPSA-Earth Authors
+#
+# SPDX-License-Identifier: AGPL-3.0-or-later
+# -*- coding: utf-8 -*-
+"""
+This script processes global building data to calculate solar rooftop area.
+"""
+import country_converter as coco
+import geopandas as gpd
+import numpy as np
+import pandas as pd
+from _helpers import (
+    BASE_DIR,
+    configure_logging,
+    create_logger,
+)
+from shapely import Point
+
+cc = coco.CountryConverter()
+
+logger = create_logger(__name__)
+
+
+def calculate_solar_rooftop_area(
+    df,
+    country_code,
+    shapes,
+    output,
+    crs,
+    install_ratio,
+    tolerance=100,
+):
+    """
+    Calculates the usable solar rooftop area for buildings within a country, considering
+    that only a portion of each rooftop can be allocated for PV installation based on its size.
+    Smaller buildings are less likely to host PV systems, while larger buildings can dedicate
+    a greater share of their rooftop to PVs. This approach is based on the methodology
+    described in a paper by Hideaki Obane (https://eneken.ieej.or.jp/data/12710.pdf).
+
+    Parameters
+    ----------
+    df : pandas.DataFrame
+        DataFrame containing building areas and center points.
+    country_code : str
+        The ISO2 country code.
+    shapes : geopandas.GeoDataFrame
+        GeoDataFrame of regional shapes with 'country' and 'geometry' columns.
+    output : str
+        The file path to save the calculated solar rooftop area.
+    crs : dict
+        A dictionary containing coordinate reference systems (CRS) for 'distance_crs' and 'geo_crs'.
+    install_ratio : dict
+        A dictionary mapping building area thresholds to solar installation ratios.
+    tolerance : int, optional
+        Maximum distance in kilometers for spatial join to find nearest shapes.
+        The default is 100.
+    """
+    distance_crs = crs["distance_crs"]
+    geo_crs = crs["geo_crs"]
+
+    keys = np.array(sorted(install_ratio.keys()))
+    values = np.array([install_ratio[k] for k in keys])
+
+    def get_ratio(area):
+        # find the largest key <= area
+        idx = np.searchsorted(keys, area, side="right") - 1
+        if idx < 0:
+            return np.nan  # or 0 if you prefer a default
+        return values[idx]
+
+    df["usefull_area"] = df["area"] * df["area"].apply(get_ratio)
+    gdf = gpd.GeoDataFrame(df, geometry="center", crs=geo_crs).to_crs(distance_crs)
+    shapes_country = shapes[shapes.country == country_code].to_crs(distance_crs)
+
+    joined = gpd.sjoin(gdf, shapes_country, how="left", predicate="intersects")
+    unmatched = joined[joined.name.isna()]
+    unmatched = unmatched.drop(["name", "country"], axis=1)
+
+    if not unmatched.empty:
+        nearest = gpd.sjoin_nearest(
+            unmatched,
+            shapes_country,
+            max_distance=tolerance * 1e3,
+        )
+
+        # Replace the unmatched rows in `joined` with the nearest results
+        matched = joined[joined.name.notna()]
+        joined = pd.concat([matched, nearest], ignore_index=True)
+
+    usefull_area = joined.groupby("name")["usefull_area"].sum()
+    shapes_country["usefull_area"] = shapes_country.index.map(usefull_area)
+
+    # Save file
+    shapes_country["usefull_area"].to_csv(output)
+
+
+if __name__ == "__main__":
+    if "snakemake" not in globals():
+        from _helpers import mock_snakemake
+
+        snakemake = mock_snakemake(
+            "cluster_global_buildings", simpl="", clusters=10, country="NG"
+        )
+
+    configure_logging(snakemake)
+    country_code = snakemake.wildcards.country
+
+    # Retrieve files
+    logger.info(f"Reading Global Buildings for {country_code}")
+    df = pd.read_parquet(snakemake.input.country_buildings)
+    df["center"] = df.apply(lambda row: Point(row["x"], row["y"]), axis=1)
+
+    shapes = gpd.read_file(snakemake.input.regions_onshore).set_index("name")[
+        ["country", "geometry"]
+    ]
+
+    if snakemake.params.solar_rooftop_enable:
+        logger.info(f"Calculate solar rooftop area for {country_code}")
+
+        calculate_solar_rooftop_area(
+            df,
+            country_code,
+            shapes,
+            snakemake.output.solar_rooftop_layout,
+            crs=snakemake.params.crs,
+            install_ratio=snakemake.params.install_ratio,
+            tolerance=snakemake.params.tolerance,
+        )