Build renewable potential

Functions associated with the build_renewable_potential rule. - Temporal Profiles are built based on the atlite cutout - Potentials are built based on the atlite cutout and raster data (land availability)

make_offshore_wind_profile(offwind_config, cutout, outp_path)

Make the offwind geographical potentials and per unit availability time series for each raster cell ! Somewhat compute intensive !

Parameters:

Name	Type	Description	Default
offwind_config	dict	the configuration for the offshore wind	required
cutout	Cutout	the atlite cutout	required
outp_path	PathLike	the output path for the raster date	required

Source code in workflow/scripts/build_renewable_potential.py

def make_offshore_wind_profile(offwind_config: dict, cutout: atlite.Cutout, outp_path: PathLike):
    """Make the offwind geographical potentials and per unit availability time series for
      each raster cell
    ! Somewhat compute intensive !


    Args:
        offwind_config (dict): the configuration for the offshore wind
        cutout (atlite.Cutout): the atlite cutout
        outp_path (PathLike): the output path for the raster date
    """
    offwind_resource = offwind_config["resource"]
    offwind_correction_factor = offwind_config.get(
        "correction_factor", DEFAULT_OFFSHORE_WIND_CORR_FACTOR
    )
    offwind_capacity_per_sqkm = offwind_config["capacity_per_sqkm"]
    if offwind_correction_factor != 1.0:
        logger.info(f"offwind_correction_factor is set as {offwind_correction_factor}")

    offwind_provinces = OFFSHORE_WIND_NODES

    EEZ_province_shp = gpd.read_file(snakemake.input["offshore_province_shapes"]).set_index(
        "province"
    )
    EEZ_province_shp = EEZ_province_shp.reindex(offwind_provinces).rename_axis("bus")
    if EEZ_province_shp.geometry.isnull().any():
        empty_geoms = EEZ_province_shp[EEZ_province_shp.geometry.isnull()].index.to_list()
        raise ValueError(
            f"There are empty geometries in EEZ_province_shp {empty_geoms}, offshore wind will fail"
        )
    EEZ_country = gpd.GeoDataFrame(
        geometry=[EEZ_province_shp.unary_union],
        crs=EEZ_province_shp.crs,
        index=["country"],
    )

    excluder_offwind = ExclusionContainer(crs=3035, res=500)

    if "max_depth" in offwind_config:
        func = functools.partial(np.greater, -offwind_config["max_depth"])
        excluder_offwind.add_raster(snakemake.input.gebco, codes=func, crs=CRS, nodata=-1000)

    if offwind_config["natura"]:
        nat1 = gpd.read_file(snakemake.input["natura1"])
        nat2 = gpd.read_file(snakemake.input["natura2"])
        nat3 = gpd.read_file(snakemake.input["natura3"])

        protected_shp = gpd.GeoDataFrame(pd.concat([nat1, nat2, nat3], ignore_index=True))
        protected_shp = gpd.GeoDataFrame(protected_shp.geometry)

        protected_Marine_shp = gpd.tools.overlay(protected_shp, EEZ_country, how="intersection")

        # this is to avoid atlite complaining about parallelisation
        logger.info("Creating tmp directory for protected marine shapefile")
        logger.info(f"parent exists: {os.path.isdir(os.path.dirname(os.path.dirname(TMP)))}")
        if not os.path.isdir(os.path.dirname(os.path.dirname(TMP))):
            mkdir(os.path.dirname(os.path.dirname(TMP)))
        if not os.path.isdir(os.path.dirname(TMP)):
            mkdir(os.path.dirname(TMP))
        protected_Marine_shp.to_file(TMP)
        excluder_offwind.add_geometry(TMP)

    kwargs = dict(nprocesses=nprocesses, disable_progressbar=noprogress)
    if noprogress:
        logger.info("Calculate offwind landuse availabilities...")
        start = time.time()
        offwind_matrix = cutout.availabilitymatrix(EEZ_province_shp, excluder_offwind, **kwargs)
        duration = time.time() - start
        logger.info(f"Completed offwind availability calculation ({duration:2.2f}s)")
    else:
        offwind_matrix = cutout.availabilitymatrix(EEZ_province_shp, excluder_offwind, **kwargs)

    offwind_potential = offwind_capacity_per_sqkm * offwind_matrix.sum("bus") * area

    offwind_func = getattr(cutout, offwind_resource.pop("method"))
    offwind_resource["dask_kwargs"] = {"num_workers": nprocesses}  # ?
    offwind_capacity_factor = offwind_correction_factor * offwind_func(
        capacity_factor=True, **offwind_resource
    )
    offwind_layout = offwind_capacity_factor * area * offwind_capacity_per_sqkm
    offwind_profile, offwind_capacities = offwind_func(
        matrix=offwind_matrix.stack(spatial=["y", "x"]),
        layout=offwind_layout,
        index=EEZ_province_shp.index,
        per_unit=True,
        return_capacity=True,
        **offwind_resource,
    )

    logger.info("Calculating offwind maximal capacity per bus (method 'simple')")

    offwind_p_nom_max = offwind_capacity_per_sqkm * offwind_matrix @ area

    offwind_ds = xr.merge(
        [
            (offwind_correction_factor * offwind_profile).rename("profile"),
            offwind_capacities.rename("weight"),
            offwind_p_nom_max.rename("p_nom_max"),
            offwind_potential.rename("potential"),
        ]
    )

    offwind_ds = offwind_ds.sel(
        bus=(
            (offwind_ds["profile"].mean("time") > offwind_config.get("min_p_max_pu", 0.0))
            & (offwind_ds["p_nom_max"] > offwind_config.get("min_p_nom_max", 0.0))
        )
    )

    if "clip_p_max_pu" in offwind_config:
        min_p_max_pu = offwind_config["clip_p_max_pu"]
        offwind_ds["profile"] = offwind_ds["profile"].where(
            offwind_ds["profile"] >= min_p_max_pu, 0
        )
    # shift back from UTC to network time
    offwind_ds["time"] = (
        pd.DatetimeIndex(offwind_ds["time"], tz="UTC").tz_convert(TIMEZONE).tz_localize(None).values
    )
    offwind_ds.to_netcdf(outp_path)

make_onshore_wind_profile(onwind_config, cutout, outp_path)

Make the onwind geographical potentials and per unit availability time series for each raster cell ! Somewhat compute intensive !

Parameters:

Name	Type	Description	Default
onwind_config	dict	the onshore wind config (from the yaml config read by snakemake)	required
cutout	Cutout	the atlite cutout	required
outp_path	PathLike	the output path for the raster data	required

Source code in workflow/scripts/build_renewable_potential.py

def make_onshore_wind_profile(onwind_config: dict, cutout: atlite.Cutout, outp_path: PathLike):
    """Make the onwind geographical potentials and per unit availability time series for
    each raster cell
    ! Somewhat compute intensive !

    Args:
        onwind_config (dict): the onshore wind config (from the yaml config read by snakemake)
        cutout (atlite.Cutout): the atlite cutout
        outp_path (PathLike): the output path for the raster data
    """

    logger.info("Making onshore wind profile ")

    onwind_resource = onwind_config["resource"]
    onwind_correction_factor = onwind_config.get("correction_factor", 1.0)
    onwind_capacity_per_sqkm = onwind_config["capacity_per_sqkm"]
    if onwind_correction_factor != 1.0:
        logger.info(f"onwind_correction_factor is set as {onwind_correction_factor}")

    excluder_onwind = ExclusionContainer(crs=3035, res=500)

    excluder_onwind.add_raster(grass, invert=True, crs=4326)
    excluder_onwind.add_raster(bare, invert=True, crs=4326)
    excluder_onwind.add_raster(shrubland, invert=True, crs=4326)

    kwargs = dict(nprocesses=nprocesses, disable_progressbar=noprogress)
    if noprogress:
        logger.info("Calculate onwind landuse availabilities...")
        start = time.time()
        onwind_matrix = cutout.availabilitymatrix(provinces_shp, excluder_onwind, **kwargs)
        duration = time.time() - start
        logger.info(f"Completed onwind availability calculation ({duration:2.2f}s)")
    else:
        onwind_matrix = cutout.availabilitymatrix(provinces_shp, excluder_onwind, **kwargs)

    onwind_potential = onwind_capacity_per_sqkm * onwind_matrix.sum("bus") * area

    onwind_func = getattr(cutout, onwind_resource.pop("method"))
    onwind_resource["dask_kwargs"] = {"num_workers": nprocesses}  # ?
    onwind_capacity_factor = onwind_correction_factor * onwind_func(
        capacity_factor=True, **onwind_resource
    )
    onwind_layout = onwind_capacity_factor * area * onwind_capacity_per_sqkm
    onwind_profile, onwind_capacities = onwind_func(
        matrix=onwind_matrix.stack(spatial=["y", "x"]),
        layout=onwind_layout,
        index=buses,
        per_unit=True,
        return_capacity=True,
        **onwind_resource,
    )

    logger.info("Calculating onwind maximal capacity per bus (method 'simple')")

    onwind_p_nom_max = onwind_capacity_per_sqkm * onwind_matrix @ area

    onwind_ds = xr.merge(
        [
            (onwind_correction_factor * onwind_profile).rename("profile"),
            onwind_capacities.rename("weight"),
            onwind_p_nom_max.rename("p_nom_max"),
            onwind_potential.rename("potential"),
        ]
    )

    onwind_ds = onwind_ds.sel(
        bus=(
            (onwind_ds["profile"].mean("time") > onwind_config.get("min_p_max_pu", 0.0))
            & (onwind_ds["p_nom_max"] > onwind_config.get("min_p_nom_max", 0.0))
        )
    )

    if "clip_p_max_pu" in onwind_config:
        min_p_max_pu = onwind_config["clip_p_max_pu"]
        onwind_ds["profile"] = onwind_ds["profile"].where(onwind_ds["profile"] >= min_p_max_pu, 0)

    # shift back from UTC to network time
    onwind_ds["time"] = (
        pd.DatetimeIndex(onwind_ds["time"], tz="UTC").tz_convert(TIMEZONE).tz_localize(None).values
    )
    onwind_ds.to_netcdf(outp_path)

make_solar_profile(solar_config, cutout, outp_path)

Make the solar geographical potentials and per unit availability time series for each raster cell ! Somewhat compute intensive !

Parameters:

Name	Type	Description	Default
solar_config	dict	the solar configuration (from the yaml config read by snakemake)	required
cutout	Cutout	the atlite cutout	required
outp_path	PathLike	the output path for the raster data	required

Source code in workflow/scripts/build_renewable_potential.py

def make_solar_profile(
    solar_config: dict,
    cutout: atlite.Cutout,
    outp_path: PathLike,
):
    """Make the solar geographical potentials and per unit availability time series for each
    raster cell
    ! Somewhat compute intensive !

    Args:
        solar_config (dict): the solar configuration (from the yaml config read by snakemake)
        cutout (atlite.Cutout): the atlite cutout
        outp_path (PathLike): the output path for the raster data
    """

    logger.info("Making solar profile ")
    solar_config = snakemake.config["renewable"]["solar"]
    solar_resource = solar_config["resource"]
    solar_correction_factor = solar_config.get("correction_factor", 1.0)
    solar_capacity_per_sqkm = solar_config["capacity_per_sqkm"]
    if solar_correction_factor != 1.0:
        logger.info(f"solar_correction_factor is set as {solar_correction_factor}")

    # TODO not hardcoded res
    excluder_solar = ExclusionContainer(crs=3035, res=500)
    excluder_build_up = ExclusionContainer(crs=3035, res=500)

    build_up = snakemake.input["Build_up_raster"]

    excluder_build_up.add_raster(build_up, invert=True, crs=CRS)
    excluder_solar.add_raster(grass, invert=True, crs=CRS)
    excluder_solar.add_raster(bare, invert=True, crs=CRS)
    excluder_solar.add_raster(shrubland, invert=True, crs=CRS)

    kwargs = dict(nprocesses=nprocesses, disable_progressbar=noprogress)
    # TODO remove if else?
    if noprogress:
        logger.info("Calculate solar landuse availabilities...")
        start = time.time()
        solar_matrix = cutout.availabilitymatrix(provinces_shp, excluder_solar, **kwargs)
        buildup_matrix = cutout.availabilitymatrix(provinces_shp, excluder_build_up, **kwargs)
        duration = time.time() - start
        logger.info(f"Completed solar availability calculation ({duration:2.2f}s)")
    else:
        solar_matrix = cutout.availabilitymatrix(
            shapes=provinces_shp, excluder=excluder_solar, **kwargs
        )
        buildup_matrix = cutout.availabilitymatrix(provinces_shp, excluder_build_up, **kwargs)

    solar_potential = (
        solar_capacity_per_sqkm * solar_matrix.sum("bus") * area
        + solar_capacity_per_sqkm * buildup_matrix.sum("bus") * area
    )

    solar_func = getattr(cutout, solar_resource.pop("method"))
    solar_resource["dask_kwargs"] = {"num_workers": nprocesses}  # ?
    solar_capacity_factor = solar_correction_factor * solar_func(
        capacity_factor=True, **solar_resource
    )
    solar_layout = solar_capacity_factor * area * solar_capacity_per_sqkm
    solar_profile, solar_capacities = solar_func(
        matrix=solar_matrix.stack(spatial=["y", "x"]),
        layout=solar_layout,
        index=buses,
        per_unit=True,
        return_capacity=True,
        **solar_resource,
    )

    logger.info("Calculating solar maximal capacity per bus (method 'simple')")

    solar_p_nom_max = solar_capacity_per_sqkm * solar_matrix @ area

    solar_ds = xr.merge(
        [
            (solar_correction_factor * solar_profile).rename("profile"),
            solar_capacities.rename("weight"),
            solar_p_nom_max.rename("p_nom_max"),
            solar_potential.rename("potential"),
        ]
    )

    solar_ds = solar_ds.sel(
        bus=(
            (solar_ds["profile"].mean("time") > solar_config.get("min_p_max_pu", 0.0))
            & (solar_ds["p_nom_max"] > solar_config.get("min_p_nom_max", 0.0))
        )
    )

    if "clip_p_max_pu" in solar_config:
        min_p_max_pu = solar_config["clip_p_max_pu"]
        solar_ds["profile"] = solar_ds["profile"].where(solar_ds["profile"] >= min_p_max_pu, 0)

    # shift back from UTC to network time
    solar_ds["time"] = (
        pd.DatetimeIndex(solar_ds["time"], tz="UTC").tz_convert(TIMEZONE).tz_localize(None).values
    )

    solar_ds.to_netcdf(outp_path)