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, delta_t)
Make the offwind geographical potentials and per unit availability time series for each raster cell ! Somewhat compute intensive !
Parameters: |
|
---|
Source code in workflow/scripts/build_renewable_potential.py
def make_offshore_wind_profile(
offwind_config: dict, cutout: atlite.Cutout, outp_path: PathLike, delta_t: pd.Timedelta
):
"""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
delta_t (pd.Timedelta): the time delta to convert to ntwk time
"""
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_shp = gpd.read_file(snakemake.input["offshore_shapes"])
EEZ_province_shp = gpd.read_file(snakemake.input["offshore_province_shapes"]).set_index("index")
EEZ_province_shp = EEZ_province_shp.reindex(offwind_provinces).rename_axis("bus")
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"]:
Protected_shp = gpd.read_file(snakemake.input["natura1"])
Protected_shp1 = gpd.read_file(snakemake.input["natura2"])
Protected_shp2 = gpd.read_file(snakemake.input["natura3"])
Protected_shp = pd.concat([Protected_shp, Protected_shp1], ignore_index=True)
Protected_shp = pd.concat([Protected_shp, Protected_shp2], ignore_index=True)
Protected_shp = Protected_shp.geometry
Protected_shp = gpd.GeoDataFrame(Protected_shp)
Protected_Marine_shp = gpd.tools.overlay(Protected_shp, EEZ_shp, how="intersection")
# this is to avoid atlite complaining about parallelisation
Protected_Marine_shp.to_file(snakemake.output.protected_areas_offshore)
# excluder_offwind.add_geometry(Protected_Marine_shp.geometry)
excluder_offwind.add_geometry("Protected_Marine.shp")
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"] = offwind_ds["time"].values + delta_t
offwind_ds.to_netcdf(outp_path)
make_onshore_wind_profile(onwind_config, cutout, outp_path, delta_t)
Make the onwind geographical potentials and per unit availability time series for each raster cell ! Somewhat compute intensive !
Parameters: |
|
---|
Source code in workflow/scripts/build_renewable_potential.py
def make_onshore_wind_profile(
onwind_config: dict, cutout: atlite.Cutout, outp_path: PathLike, delta_t: pd.Timedelta
):
"""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
delta_t (pd.Timedelta): the time delta to convert to ntwk time
"""
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"] = onwind_ds["time"].values + delta_t
onwind_ds.to_netcdf(outp_path)
make_solar_profile(solar_config, cutout, outp_path, delta_t)
Make the solar geographical potentials and per unit availability time series for each raster cell ! Somewhat compute intensive !
Parameters: |
|
---|
Source code in workflow/scripts/build_renewable_potential.py
def make_solar_profile(
solar_config: dict,
cutout: atlite.Cutout,
outp_path: PathLike,
delta_t: pd.Timedelta,
):
"""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
delta_t (pd.Timedelta): the time delta to convert to ntwk time
"""
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"] = solar_ds["time"].values + delta_t
solar_ds.to_netcdf(outp_path)