Use default cutout for ship and nature raster extents

This avoids the need to either collect all cutouts for all weather
years or arbitrarily pick all cutouts for a certain weather year.
This commit is contained in:
Koen van Greevenbroek 2024-03-01 11:45:39 +01:00
parent de0c0cd1f4
commit f6f98c102c
3 changed files with 12 additions and 22 deletions

View File

@ -205,26 +205,15 @@ if config["enable"].get("build_cutout", False):
"../scripts/build_cutout.py" "../scripts/build_cutout.py"
# TODO: consider if it's fine to just use the default cutout instead, i.e.:
# cutout="cutouts/" + CDIR + config["atlite"]["default_cutout"] + ".nc",
def cutouts_for_extent_input(wildcards):
# We only need these cutouts in order to determine geographic
# boundaries, so just pick the first weather year available
weather_year = config_provider("scenario", "weather_year")(wildcards)[0]
cutouts = []
for k in config_provider("electricity", "renewable_carriers")(wildcards):
cutout = config_provider("renewable", k, "cutout")(wildcards)
cutout = cutout.replace("{weather_year}", weather_year)
cutouts.append(f"cutouts/{CDIR}{cutout}.nc")
return cutouts
if config["enable"].get("build_natura_raster", False): if config["enable"].get("build_natura_raster", False):
rule build_natura_raster: rule build_natura_raster:
input: input:
natura=ancient("data/bundle/natura/Natura2000_end2015.shp"), natura=ancient("data/bundle/natura/Natura2000_end2015.shp"),
cutouts=cutouts_for_extent_input, cutout=lambda w: "cutouts/"
+ CDIR
+ config_provider("atlite", "default_cutout")(w)
+ ".nc",
output: output:
resources("natura.tiff"), resources("natura.tiff"),
resources: resources:
@ -240,7 +229,10 @@ if config["enable"].get("build_natura_raster", False):
rule build_ship_raster: rule build_ship_raster:
input: input:
ship_density="data/shipdensity_global.zip", ship_density="data/shipdensity_global.zip",
cutouts=cutouts_for_extent_input, cutout=lambda w: "cutouts/"
+ CDIR
+ config_provider("atlite", "default_cutout")(w)
+ ".nc",
output: output:
resources("shipdensity_raster.tif"), resources("shipdensity_raster.tif"),
log: log:

View File

@ -94,9 +94,8 @@ if __name__ == "__main__":
configure_logging(snakemake) configure_logging(snakemake)
set_scenario_config(snakemake) set_scenario_config(snakemake)
cutouts = snakemake.input.cutouts x, X, y, Y = determine_cutout_xXyY(snakemake.input.cutout)
xs, Xs, ys, Ys = zip(*(determine_cutout_xXyY(cutout) for cutout in cutouts)) bounds = transform_bounds(4326, 3035, x, y, X, Y)
bounds = transform_bounds(4326, 3035, min(xs), min(ys), max(Xs), max(Ys))
transform, out_shape = get_transform_and_shape(bounds, res=100) transform, out_shape = get_transform_and_shape(bounds, res=100)
# adjusted boundaries # adjusted boundaries

View File

@ -59,8 +59,7 @@ if __name__ == "__main__":
configure_logging(snakemake) configure_logging(snakemake)
set_scenario_config(snakemake) set_scenario_config(snakemake)
cutouts = snakemake.input.cutouts x, X, y, Y = determine_cutout_xXyY(snakemake.input.cutout)
xs, Xs, ys, Ys = zip(*(determine_cutout_xXyY(cutout) for cutout in cutouts))
with zipfile.ZipFile(snakemake.input.ship_density) as zip_f: with zipfile.ZipFile(snakemake.input.ship_density) as zip_f:
resources = Path(snakemake.output[0]).parent resources = Path(snakemake.output[0]).parent
@ -68,7 +67,7 @@ if __name__ == "__main__":
zip_f.extract(fn, resources) zip_f.extract(fn, resources)
with rioxarray.open_rasterio(resources / fn) as ship_density: with rioxarray.open_rasterio(resources / fn) as ship_density:
ship_density = ship_density.drop_vars(["band"]).sel( ship_density = ship_density.drop_vars(["band"]).sel(
x=slice(min(xs), max(Xs)), y=slice(max(Ys), min(ys)) x=slice(x, X), y=slice(Y, y)
) )
ship_density.rio.to_raster(snakemake.output[0]) ship_density.rio.to_raster(snakemake.output[0])