diff --git a/workflow/scripts/intersect/plot_wind_fields.py b/src/open_gira/wind_plotting.py similarity index 100% rename from workflow/scripts/intersect/plot_wind_fields.py rename to src/open_gira/wind_plotting.py diff --git a/tests/integration/estimate_wind_fields/data/results/power/by_country/PRI/storms/IBTrACS/tracks.geoparquet b/tests/integration/estimate_wind_fields/data/results/power/by_country/PRI/storms/IBTrACS/0/tracks.geoparquet similarity index 100% rename from tests/integration/estimate_wind_fields/data/results/power/by_country/PRI/storms/IBTrACS/tracks.geoparquet rename to tests/integration/estimate_wind_fields/data/results/power/by_country/PRI/storms/IBTrACS/0/tracks.geoparquet diff --git a/tests/integration/estimate_wind_fields/data/results/power/by_country/PRI/storms/downscale_factors.npy b/tests/integration/estimate_wind_fields/data/results/power/by_country/PRI/storms/downscale_factors.npy new file mode 100644 index 00000000..03213b64 Binary files /dev/null and b/tests/integration/estimate_wind_fields/data/results/power/by_country/PRI/storms/downscale_factors.npy differ diff --git a/tests/integration/estimate_wind_fields/data/results/power/by_country/PRI/storms/surface_roughness.tiff b/tests/integration/estimate_wind_fields/data/results/power/by_country/PRI/storms/surface_roughness.tiff deleted file mode 100644 index 975e8203..00000000 Binary files a/tests/integration/estimate_wind_fields/data/results/power/by_country/PRI/storms/surface_roughness.tiff and /dev/null differ diff --git a/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/max_wind_field.nc b/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/max_wind_field.nc similarity index 100% rename from tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/max_wind_field.nc rename to tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/max_wind_field.nc diff --git a/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/plots/2001235N18296_max_contour.png b/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2001235N18296_max_contour.png similarity index 100% rename from tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/plots/2001235N18296_max_contour.png rename to tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2001235N18296_max_contour.png diff --git a/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/plots/2006213N16302_max_contour.png b/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2006213N16302_max_contour.png similarity index 100% rename from tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/plots/2006213N16302_max_contour.png rename to tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2006213N16302_max_contour.png diff --git a/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/plots/2007345N18298_max_contour.png b/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2007345N18298_max_contour.png similarity index 100% rename from tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/plots/2007345N18298_max_contour.png rename to tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2007345N18298_max_contour.png diff --git a/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/plots/2008229N18293_max_contour.png b/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2008229N18293_max_contour.png similarity index 100% rename from tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/plots/2008229N18293_max_contour.png rename to tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2008229N18293_max_contour.png diff --git a/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/plots/2008287N15291_max_contour.png b/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2008287N15291_max_contour.png similarity index 100% rename from tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/plots/2008287N15291_max_contour.png rename to tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2008287N15291_max_contour.png diff --git a/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/plots/2009245N17303_max_contour.png b/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2009245N17303_max_contour.png similarity index 100% rename from tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/plots/2009245N17303_max_contour.png rename to tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2009245N17303_max_contour.png diff --git a/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/plots/2010236N12341_max_contour.png b/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2010236N12341_max_contour.png similarity index 100% rename from tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/plots/2010236N12341_max_contour.png rename to tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2010236N12341_max_contour.png diff --git a/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/plots/2010244N12328_max_contour.png b/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2010244N12328_max_contour.png similarity index 100% rename from tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/plots/2010244N12328_max_contour.png rename to tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2010244N12328_max_contour.png diff --git a/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/plots/2011214N15299_max_contour.png b/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2011214N15299_max_contour.png similarity index 100% rename from tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/plots/2011214N15299_max_contour.png rename to tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2011214N15299_max_contour.png diff --git a/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/plots/2011233N15301_max_contour.png b/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2011233N15301_max_contour.png similarity index 100% rename from tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/plots/2011233N15301_max_contour.png rename to tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2011233N15301_max_contour.png diff --git a/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/plots/2011250N12324_max_contour.png b/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2011250N12324_max_contour.png similarity index 100% rename from tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/plots/2011250N12324_max_contour.png rename to tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2011250N12324_max_contour.png diff --git a/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/plots/2012287N15297_max_contour.png b/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2012287N15297_max_contour.png similarity index 100% rename from tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/plots/2012287N15297_max_contour.png rename to tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2012287N15297_max_contour.png diff --git a/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/plots/2013248N16294_max_contour.png b/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2013248N16294_max_contour.png similarity index 100% rename from tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/plots/2013248N16294_max_contour.png rename to tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2013248N16294_max_contour.png diff --git a/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/plots/2014210N10323_max_contour.png b/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2014210N10323_max_contour.png similarity index 100% rename from tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/plots/2014210N10323_max_contour.png rename to tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2014210N10323_max_contour.png diff --git a/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/plots/2014285N16305_max_contour.png b/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2014285N16305_max_contour.png similarity index 100% rename from tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/plots/2014285N16305_max_contour.png rename to tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2014285N16305_max_contour.png diff --git a/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/plots/2015237N14315_max_contour.png b/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2015237N14315_max_contour.png similarity index 100% rename from tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/plots/2015237N14315_max_contour.png rename to tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2015237N14315_max_contour.png diff --git a/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/plots/2017242N16333_max_contour.png b/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2017242N16333_max_contour.png similarity index 100% rename from tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/plots/2017242N16333_max_contour.png rename to tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2017242N16333_max_contour.png diff --git a/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/plots/2017260N12310_max_contour.png b/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2017260N12310_max_contour.png similarity index 100% rename from tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/plots/2017260N12310_max_contour.png rename to tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2017260N12310_max_contour.png diff --git a/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/plots/2018186N10325_max_contour.png b/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2018186N10325_max_contour.png similarity index 100% rename from tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/plots/2018186N10325_max_contour.png rename to tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2018186N10325_max_contour.png diff --git a/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/plots/2019236N10314_max_contour.png b/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2019236N10314_max_contour.png similarity index 100% rename from tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/plots/2019236N10314_max_contour.png rename to tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2019236N10314_max_contour.png diff --git a/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/plots/2019265N12301_max_contour.png b/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2019265N12301_max_contour.png similarity index 100% rename from tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/plots/2019265N12301_max_contour.png rename to tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2019265N12301_max_contour.png diff --git a/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/plots/2020211N13306_max_contour.png b/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2020211N13306_max_contour.png similarity index 100% rename from tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/plots/2020211N13306_max_contour.png rename to tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2020211N13306_max_contour.png diff --git a/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/plots/2020233N14313_max_contour.png b/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2020233N14313_max_contour.png similarity index 100% rename from tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/plots/2020233N14313_max_contour.png rename to tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2020233N14313_max_contour.png diff --git a/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/plots/2021222N14301_max_contour.png b/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2021222N14301_max_contour.png similarity index 100% rename from tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/plots/2021222N14301_max_contour.png rename to tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2021222N14301_max_contour.png diff --git a/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/plots/2021225N15313_max_contour.png b/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2021225N15313_max_contour.png similarity index 100% rename from tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/plots/2021225N15313_max_contour.png rename to tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2021225N15313_max_contour.png diff --git a/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/plots/2022246N18301_max_contour.png b/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2022246N18301_max_contour.png similarity index 100% rename from tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/plots/2022246N18301_max_contour.png rename to tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2022246N18301_max_contour.png diff --git a/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/plots/2022257N16312_max_contour.png b/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2022257N16312_max_contour.png similarity index 100% rename from tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/plots/2022257N16312_max_contour.png rename to tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/0/plots/2022257N16312_max_contour.png diff --git a/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/downscale_factors.png b/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/downscale_factors.png deleted file mode 100644 index caac318b..00000000 Binary files a/tests/integration/estimate_wind_fields/expected/results/power/by_country/PRI/storms/IBTrACS/downscale_factors.png and /dev/null differ diff --git a/tests/integration/test_estimate_wind_fields.py b/tests/integration/test_estimate_wind_fields.py index 6bc3c923..b8227c30 100644 --- a/tests/integration/test_estimate_wind_fields.py +++ b/tests/integration/test_estimate_wind_fields.py @@ -5,6 +5,6 @@ def test_estimate_wind_fields(): runner.run_snakemake_test( "estimate_wind_fields", ( - "results/power/by_country/PRI/storms/IBTrACS/max_wind_field.nc", + "results/power/by_country/PRI/storms/IBTrACS/0/max_wind_field.nc", ) ) diff --git a/workflow/rules/exposure/electricity_grid/intersection.smk b/workflow/rules/exposure/electricity_grid/intersection.smk index 945b084b..acfe7b9d 100644 --- a/workflow/rules/exposure/electricity_grid/intersection.smk +++ b/workflow/rules/exposure/electricity_grid/intersection.smk @@ -134,7 +134,7 @@ rule electricity_grid_damages: """ input: grid_splits = rules.rasterise_electricity_grid.output.geoparquet, - wind_speeds = rules.estimate_wind_fields.output.wind_speeds, + wind_speeds = rules.concat_wind_fields_over_sample.output.concat, grid_edges = rules.create_power_network.output.edges, grid_nodes = rules.create_power_network.output.nodes, output: diff --git a/workflow/rules/exposure/wind_fields.smk b/workflow/rules/exposure/wind_fields.smk index 10f6f230..84ba5092 100644 --- a/workflow/rules/exposure/wind_fields.smk +++ b/workflow/rules/exposure/wind_fields.smk @@ -156,11 +156,29 @@ snakemake -c1 results/power/by_country/CUB/surface_roughness.tiff """ +rule create_downscaling_factors: + """ + Use surface roughness values calculate factors for scaling from gradient-level to surface-level winds. + """ + input: + surface_roughness=rules.create_surface_roughness_raster.output.surface_roughness, + output: + downscale_factors="{OUTPUT_DIR}/power/by_country/{COUNTRY_ISO_A3}/storms/downscale_factors.npy", + downscale_factors_plot="{OUTPUT_DIR}/power/by_country/{COUNTRY_ISO_A3}/storms/downscale_factors.png", + script: + "../../scripts/intersect/wind_downscaling_factors.py" + +""" +To test: +snakemake -c1 results/power/by_country/PRI/storms/downscale_factors.npy +""" + + def storm_tracks_by_country(wildcards) -> str: """Return path to storm tracks""" # parent dataset of storm set e.g. IBTrACS for IBTrACS_maria-2017 storm_dataset = wildcards.STORM_SET.split("_")[0] - return f"{wildcards.OUTPUT_DIR}/power/by_country/{wildcards.COUNTRY_ISO_A3}/storms/{storm_dataset}/tracks.geoparquet" + return f"{wildcards.OUTPUT_DIR}/power/by_country/{wildcards.COUNTRY_ISO_A3}/storms/{storm_dataset}/{wildcards.SAMPLE}/tracks.geoparquet" def read_storm_set(wildcards) -> list[str]: @@ -185,24 +203,101 @@ rule estimate_wind_fields: input: storm_file=storm_tracks_by_country, wind_grid="{OUTPUT_DIR}/power/by_country/{COUNTRY_ISO_A3}/storms/wind_grid.tiff", - surface_roughness=rules.create_surface_roughness_raster.output.surface_roughness, + downscaling_factors=rules.create_downscaling_factors.output.downscale_factors, params: storm_set=read_storm_set, # include failure_thresholds as a param (despite not using elsewhere in the # rule) to trigger re-runs on change to this configuration option failure_thresholds=config["transmission_windspeed_failure"] threads: threads_for_country + resources: + # 2GB RAM per CPU + mem_mb = lambda wildcards: threads_for_country(wildcards) * 2_048 output: # enable or disable plotting in the config file - plot_dir=directory("{OUTPUT_DIR}/power/by_country/{COUNTRY_ISO_A3}/storms/{STORM_SET}/plots/"), - wind_speeds="{OUTPUT_DIR}/power/by_country/{COUNTRY_ISO_A3}/storms/{STORM_SET}/max_wind_field.nc", - downscale_factors_plot="{OUTPUT_DIR}/power/by_country/{COUNTRY_ISO_A3}/storms/{STORM_SET}/downscale_factors.png", + plot_dir=directory("{OUTPUT_DIR}/power/by_country/{COUNTRY_ISO_A3}/storms/{STORM_SET}/{SAMPLE}/plots/"), + wind_speeds="{OUTPUT_DIR}/power/by_country/{COUNTRY_ISO_A3}/storms/{STORM_SET}/{SAMPLE}/max_wind_field.nc", script: "../../scripts/intersect/estimate_wind_fields.py" """ To test: -snakemake -c1 results/power/by_country/PRI/storms/IBTrACS/max_wind_field.nc +snakemake -c1 results/power/by_country/PRI/storms/IBTrACS/0/max_wind_field.nc +""" + + +def wind_field_paths_all_samples(wildcards) -> list[str]: + """ + Return a list of paths for every sample + """ + dataset_name = wildcards.STORM_SET.split("-")[0] + return expand( + "{OUTPUT_DIR}/power/by_country/{COUNTRY_ISO_A3}/storms/{STORM_SET}/{SAMPLE}/max_wind_field.nc", + OUTPUT_DIR=wildcards.OUTPUT_DIR, + COUNTRY_ISO_A3=wildcards.COUNTRY_ISO_A3, + STORM_SET=wildcards.STORM_SET, + SAMPLE=range(0, SAMPLES_PER_TRACKSET[dataset_name]) + ) + + +rule concat_wind_fields_over_sample: + """ + Take wind fields expanded over sample and stack them into a single netCDF. + """ + input: + sample_paths=wind_field_paths_all_samples + output: + concat="{OUTPUT_DIR}/power/by_country/{COUNTRY_ISO_A3}/storms/{STORM_SET}/max_wind_field.nc", + run: + import logging + + import dask + import xarray as xr + + from open_gira.io import netcdf_packing_parameters + + logging.basicConfig(format="%(asctime)s %(process)d %(filename)s %(message)s", level=logging.INFO) + + logging.info("Reading wind fields from each sample") + # concatenated xarray dataset is chunked by input file + # N.B. this is lazily loaded + all_samples: xr.Dataset = xr.open_mfdataset( + input.sample_paths, + chunks={"max_wind_speed": len(input.sample_paths)}, + ) + + logging.info("Computing packing factors for all samples") + # we used dask to allow for chunked calculation of data min/max and to stream to disk + # use the synchronous scheduler to limit dask to a single process (reducing memory footprint) + scheduler = "synchronous" + + # compute packing factors for all data, need global min and max + # the implementation below reads all the data chunks twice, once for min and once for max + # would be nice to avoid this duplication + scale_factor, add_offset, fill_value = netcdf_packing_parameters( + all_samples.max_wind_speed.min().compute(scheduler=scheduler).item(), + all_samples.max_wind_speed.max().compute(scheduler=scheduler).item(), + 16 + ) + + logging.info("Writing pooled wind fields to disk") + serialisation_task: dask.delayed.Delayed = all_samples.to_netcdf( + output.concat, + encoding={ + "max_wind_speed": { + 'dtype': 'int16', + 'scale_factor': scale_factor, + 'add_offset': add_offset, + '_FillValue': fill_value + } + }, + compute=False + ) + serialisation_task.compute(scheduler=scheduler) + +""" +To test: +snakemake -c1 results/power/by_country/PRI/storms/STORM-constant/max_wind_field.nc """ @@ -216,9 +311,10 @@ def wind_fields_by_country_for_storm(wildcards): country_set_by_storm = cached_json_file_read(json_file) return expand( - "results/power/by_country/{COUNTRY_ISO_A3}/storms/{STORM_SET}/max_wind_field.nc", + "results/power/by_country/{COUNTRY_ISO_A3}/storms/{STORM_SET}/{SAMPLE}/max_wind_field.nc", COUNTRY_ISO_A3=country_set_by_storm[wildcards.STORM_ID], # list of str STORM_SET=wildcards.STORM_SET, # str + SAMPLE=wildcards.SAMPLE # str ) @@ -229,7 +325,7 @@ rule merge_wind_fields_of_storm: input: wind_fields = wind_fields_by_country_for_storm output: - merged = "{OUTPUT_DIR}/power/by_storm_set/{STORM_SET}/by_storm/{STORM_ID}/wind_field.nc", + merged = "{OUTPUT_DIR}/power/by_storm_set/{STORM_SET}/{SAMPLE}/by_storm/{STORM_ID}/wind_field.nc", run: import logging import os @@ -298,7 +394,7 @@ rule merge_wind_fields_of_storm: """ Test with: -snakemake -c1 results/power/by_storm_set/IBTrACS/by_storm/2017260N12310/wind_field.nc +snakemake -c1 results/power/by_storm_set/IBTrACS/0/by_storm/2017260N12310/wind_field.nc """ @@ -315,9 +411,10 @@ def merged_wind_fields_for_all_storms_in_storm_set(wildcards): storms = list(country_set_by_storm.keys()) return expand( - "results/power/by_storm_set/{STORM_SET}/by_storm/{STORM_ID}/wind_field.nc", + "results/power/by_storm_set/{STORM_SET}/{SAMPLE}/by_storm/{STORM_ID}/wind_field.nc", STORM_SET=wildcards.STORM_SET, # str - STORM_ID=storms # list of str + STORM_ID=storms, # list of str + SAMPLE=wildcards.SAMPLE # str ) @@ -329,7 +426,7 @@ rule merged_wind_fields_for_storm_set: input: wind_fields = merged_wind_fields_for_all_storms_in_storm_set output: - completion_flag = "{OUTPUT_DIR}/power/by_storm_set/{STORM_SET}/wind_fields.txt" + completion_flag = "{OUTPUT_DIR}/power/by_storm_set/{STORM_SET}/{SAMPLE}/wind_fields.txt" shell: """ # one output file per line @@ -338,5 +435,5 @@ rule merged_wind_fields_for_storm_set: """ Test with: -snakemake -c1 -- results/power/by_storm_set/IBTrACS/wind_fields.txt +snakemake -c1 -- results/power/by_storm_set/IBTrACS/0/wind_fields.txt """ diff --git a/workflow/rules/preprocess/IBTrACS.smk b/workflow/rules/preprocess/IBTrACS.smk index 4c17d65c..585580e5 100644 --- a/workflow/rules/preprocess/IBTrACS.smk +++ b/workflow/rules/preprocess/IBTrACS.smk @@ -2,13 +2,13 @@ rule parse_ibtracs: input: ibtracs_csv = "results/input/IBTrACS/raw/v4.csv" output: - ibtracs_parquet = "results/storm_tracks/IBTrACS/tracks.geoparquet" + ibtracs_parquet = "results/storm_tracks/IBTrACS/0/tracks.geoparquet" script: "../../scripts/preprocess/parse_IBTrACS.py" """ To test: -snakemake -c1 results/storm_tracks/IBTrACS/tracks.geoparquet +snakemake -c1 results/storm_tracks/IBTrACS/0/tracks.geoparquet """ @@ -17,11 +17,11 @@ rule slice_ibtracs: global_tracks=rules.parse_ibtracs.output.ibtracs_parquet, grid_hull="{OUTPUT_DIR}/power/by_country/{COUNTRY_ISO_A3}/network/convex_hull.json" output: - sliced_tracks="{OUTPUT_DIR}/power/by_country/{COUNTRY_ISO_A3}/storms/IBTrACS/tracks.geoparquet", + sliced_tracks="{OUTPUT_DIR}/power/by_country/{COUNTRY_ISO_A3}/storms/IBTrACS/0/tracks.geoparquet", script: "../../scripts/preprocess/slice_storm_tracks.py" """ To test: -snakemake -c1 results/power/by_country/PRI/storms/IBTrACS/tracks.geoparquet +snakemake -c1 results/power/by_country/PRI/storms/IBTrACS/0/tracks.geoparquet """ diff --git a/workflow/rules/preprocess/IRIS.smk b/workflow/rules/preprocess/IRIS.smk index 9695023c..4b69c002 100644 --- a/workflow/rules/preprocess/IRIS.smk +++ b/workflow/rules/preprocess/IRIS.smk @@ -2,13 +2,13 @@ rule parse_IRIS: input: csv_dir="{OUTPUT_DIR}/input/IRIS/iris-data/event_sets/{IRIS_SCENARIO}/" output: - parquet="{OUTPUT_DIR}/storm_tracks/IRIS-{IRIS_SCENARIO}/tracks.geoparquet" + parquet="{OUTPUT_DIR}/storm_tracks/IRIS-{IRIS_SCENARIO}/{SAMPLE}/tracks.geoparquet" script: "../../scripts/preprocess/parse_IRIS.py" """ Test with: -snakemake -c1 results/storm_tracks/IRIS_SSP1-2050/tracks.geoparquet +snakemake -c1 results/storm_tracks/IRIS_SSP1-2050/0/tracks.geoparquet """ @@ -17,7 +17,7 @@ rule slice_IRIS: global_tracks=rules.parse_IRIS.output.parquet, grid_hull="{OUTPUT_DIR}/power/by_country/{COUNTRY_ISO_A3}/network/convex_hull.json" output: - sliced_tracks="{OUTPUT_DIR}/power/by_country/{COUNTRY_ISO_A3}/storms/IRIS-{IRIS_SCENARIO}/tracks.geoparquet", + sliced_tracks="{OUTPUT_DIR}/power/by_country/{COUNTRY_ISO_A3}/storms/IRIS-{IRIS_SCENARIO}/{SAMPLE}/tracks.geoparquet", resources: mem_mb=10000 # the global tracks file is fairly chunky script: @@ -25,5 +25,5 @@ rule slice_IRIS: """ To test: -snakemake -c1 results/power/by_country/PRI/storms/IRIS-PRESENT/tracks.geoparquet +snakemake -c1 results/power/by_country/PRI/storms/IRIS-PRESENT/0/tracks.geoparquet """ diff --git a/workflow/rules/preprocess/STORM.smk b/workflow/rules/preprocess/STORM.smk index a99c9b38..abfdc79f 100644 --- a/workflow/rules/preprocess/STORM.smk +++ b/workflow/rules/preprocess/STORM.smk @@ -2,22 +2,22 @@ rule parse_storm: input: csv_dir="{OUTPUT_DIR}/input/STORM/events/{STORM_MODEL}/raw" output: - parquet="{OUTPUT_DIR}/storm_tracks/STORM-{STORM_MODEL}/tracks.geoparquet" + parquet="{OUTPUT_DIR}/storm_tracks/STORM-{STORM_MODEL}/{SAMPLE}/tracks.geoparquet" script: "../../scripts/preprocess/parse_STORM.py" """ Test with: -snakemake -c1 results/storm_tracks/STORM-constant/tracks.geoparquet +snakemake -c1 results/storm_tracks/STORM-constant/0/tracks.geoparquet """ rule slice_storm: input: - global_tracks=rules.parse_storm.output.parquet, + global_tracks="{OUTPUT_DIR}/storm_tracks/STORM-{STORM_MODEL}/{SAMPLE}/tracks.geoparquet", grid_hull="{OUTPUT_DIR}/power/by_country/{COUNTRY_ISO_A3}/network/convex_hull.json" output: - sliced_tracks="{OUTPUT_DIR}/power/by_country/{COUNTRY_ISO_A3}/storms/STORM-{STORM_MODEL}/tracks.geoparquet", + sliced_tracks="{OUTPUT_DIR}/power/by_country/{COUNTRY_ISO_A3}/storms/STORM-{STORM_MODEL}/{SAMPLE}/tracks.geoparquet", resources: mem_mb=10000 # the global tracks file is fairly chunky script: @@ -25,5 +25,5 @@ rule slice_storm: """ To test: -snakemake -c1 results/power/by_country/PRI/storms/STORM-constant/tracks.geoparquet -""" +snakemake -c1 results/power/by_country/PRI/storms/STORM-constant/0/tracks.geoparquet +""" \ No newline at end of file diff --git a/workflow/rules/preprocess/create_network.smk b/workflow/rules/preprocess/create_network.smk index e78aae1b..111c6ba0 100644 --- a/workflow/rules/preprocess/create_network.smk +++ b/workflow/rules/preprocess/create_network.smk @@ -240,7 +240,7 @@ def threads_for_country(wildcards) -> int: ranked["threads"] = logistic_min( ranked.index, # input to transform - workflow.cores, # maximum (roughly) + 0.8 * workflow.cores, # maximum (roughly) -2, # minimum (roughly) 0.02, # steepness of sigmoid 0.8 * len(ranked.index) # location of sigmoid centre on input axis diff --git a/workflow/rules/preprocess/join_data.smk b/workflow/rules/preprocess/join_data.smk index c1709097..805d383f 100644 --- a/workflow/rules/preprocess/join_data.smk +++ b/workflow/rules/preprocess/join_data.smk @@ -1,8 +1,8 @@ -# Take split .geoparquet files and output a single, unified .geoparquet file -# for each dataset-hazard combination - - rule join_data: + """ + Take split .geoparquet files and output a single, unified .geoparquet file + for each dataset-hazard combination + """ input: lambda wildcards: expand( os.path.join( @@ -20,8 +20,42 @@ rule join_data: script: "../../scripts/join_data.py" - """ Test with: snakemake --cores all results/tanzania-mini_filter-road_hazard-aqueduct-river.geoparquet """ + + +def all_storm_tracks_by_sample(wildcards) -> list[str]: + """ + Return a list of every per-sample tracks file for a given STORM_SET. + """ + dataset_name = wildcards.STORM_SET.split("-")[0] + return expand( + "{OUTPUT_DIR}/storm_tracks/{STORM_SET}/{SAMPLE}/tracks.geoparquet", + OUTPUT_DIR=wildcards.OUTPUT_DIR, + STORM_SET=wildcards.STORM_SET, + SAMPLE=range(0, SAMPLES_PER_TRACKSET[dataset_name]) + ) + + +rule concat_storm_tracks: + """ + To concat all samples of tracks together, useful when we want to e.g. plot track atlas. + """ + input: + by_sample=all_storm_tracks_by_sample + output: + tracks_from_all_samples="{OUTPUT_DIR}/storm_tracks/{STORM_SET}/tracks.geoparquet", + run: + import geopandas as gpd + import pandas as pd + + data = [gpd.read_parquet(path) for path in input.by_sample] + concat = gpd.GeoDataFrame(pd.concat(data)) + concat.to_parquet(output.tracks_from_all_samples) + +""" +Test with: +snakemake -c1 -- results/storm_tracks/STORM-constant/tracks.geoparquet +""" \ No newline at end of file diff --git a/workflow/rules/storm_workflow_global_variables.smk b/workflow/rules/storm_workflow_global_variables.smk index 386b8727..10e85446 100644 --- a/workflow/rules/storm_workflow_global_variables.smk +++ b/workflow/rules/storm_workflow_global_variables.smk @@ -32,6 +32,13 @@ def country_codes() -> List[str]: # list of ISO A3 country codes COUNTRY_CODES = country_codes() +# how many samples is each storm track dataset split into? +SAMPLES_PER_TRACKSET = { + "IBTrACS": 1, + "STORM": 10, + "IRIS": 10, +} + STORM_RPS = ( list(range(10, 100, 10)) + list(range(100, 1000, 100)) diff --git a/workflow/scripts/intersect/estimate_wind_fields.py b/workflow/scripts/intersect/estimate_wind_fields.py index 1e882901..2e8e0d03 100644 --- a/workflow/scripts/intersect/estimate_wind_fields.py +++ b/workflow/scripts/intersect/estimate_wind_fields.py @@ -19,30 +19,20 @@ from open_gira.io import netcdf_packing_parameters from open_gira.wind import ( advective_vector, rotational_field, interpolate_track, - power_law_scale_factors, empty_wind_da, WIND_COORDS, ENV_PRESSURE + empty_wind_da, WIND_COORDS, ENV_PRESSURE ) -from plot_wind_fields import plot_contours, animate_track, plot_downscale_factors +from open_gira.wind_plotting import plot_contours, animate_track logging.basicConfig(format="%(asctime)s %(process)d %(filename)s %(message)s", level=logging.INFO) -# wind speed altitudes -# gradient level clearly not realistic, but we vary it to fit our estimated wind -# speeds to observations (or to better model results) -# the 18m figure is as a result of minimising pixel-wise errors between this model -# and that used in Done et al. 2020 with a physical boundary layer -GRADIENT_LEVEL_METRES = 18 -SURFACE_LEVEL_METRES = 10 - - -def cleanup(output_path: str, downscale_factors_plot_path: str): +def cleanup(output_path: str): """ If we don't have a network, or tracks and we can't continue, write empty - files and quit. + file and quit. """ empty_wind_da().to_netcdf(output_path) - os.system(f"touch {downscale_factors_plot_path}") sys.exit(0) @@ -177,14 +167,13 @@ def process_track( storm_file_path: str = snakemake.input.storm_file wind_grid_path: str = snakemake.input.wind_grid - surface_roughness_path: str = snakemake.input.surface_roughness + downscale_factors_path: str = snakemake.input.downscaling_factors storm_set: set[str] = set(snakemake.params.storm_set) plot_max_wind: bool = snakemake.config["plot_wind"]["max_speed"] plot_animation: bool = snakemake.config["plot_wind"]["animation"] n_proc: int = snakemake.threads plot_dir_path: str = snakemake.output.plot_dir output_path: str = snakemake.output.wind_speeds - downscale_factors_plot_path: str = snakemake.output.downscale_factors_plot # directory required (snakemake output) os.makedirs(plot_dir_path, exist_ok=True) @@ -193,7 +182,7 @@ def process_track( tracks = gpd.read_parquet(storm_file_path) if tracks.empty: logging.info("No intersection between network and tracks, writing empty file.") - cleanup(output_path, downscale_factors_plot_path) + cleanup(output_path) if storm_set: # if we have a storm_set, only keep the matching track_ids @@ -202,7 +191,7 @@ def process_track( if tracks.empty: logging.info("No intersection between network and tracks, writing empty file.") - cleanup(output_path, downscale_factors_plot_path) + cleanup(output_path) logging.info(f"\n{tracks}") grouped_tracks = tracks.groupby("track_id") @@ -212,29 +201,8 @@ def process_track( grid: xr.DataArray = rioxarray.open_rasterio(wind_grid_path) logging.info(f"\n{grid}") - logging.info("Calculating downscaling factors from surface roughness raster") - - # surface roughness raster for downscaling winds with - surface_roughness_raster: xr.DataArray = rioxarray.open_rasterio(surface_roughness_path) - - # (1, y, x) where 1 is the number of surface roughness bands - # select the only value in the band dimension - surface_roughness: np.ndarray = surface_roughness_raster[0].values - - # calculate factors to scale wind speeds from gradient-level to surface level, - # taking into account surface roughness as defined by the raster - downscaling_factors = power_law_scale_factors( - surface_roughness, - SURFACE_LEVEL_METRES, - GRADIENT_LEVEL_METRES - ) - - logging.info("Mapping downscaling factors and saving to disk") - plot_downscale_factors( - downscaling_factors, - "Wind downscaling factors", - downscale_factors_plot_path - ) + logging.info("Reading wind downscaling factors") + downscaling_factors = np.load(downscale_factors_path) # track is a tuple of track_id and the tracks subset, we only want the latter args = ((track[1], grid.x, grid.y, downscaling_factors, plot_max_wind, plot_animation, plot_dir_path) for track in grouped_tracks) diff --git a/workflow/scripts/intersect/wind_downscaling_factors.py b/workflow/scripts/intersect/wind_downscaling_factors.py new file mode 100644 index 00000000..077ed271 --- /dev/null +++ b/workflow/scripts/intersect/wind_downscaling_factors.py @@ -0,0 +1,66 @@ +""" +Estimate downscaling factors from gradient to surface level winds, given surface +roughness raster. +""" + +import logging +import os + +import numpy as np +import rasterio +import rioxarray +import xarray as xr + +from open_gira.wind import power_law_scale_factors +from open_gira.wind_plotting import plot_downscale_factors + +if __name__ == "__main__": + + logging.basicConfig(format="%(asctime)s %(process)d %(filename)s %(message)s", level=logging.INFO) + + # wind speed altitudes + + # gradient level clearly not realistic, but we vary it to fit our estimated wind + # speeds to observations (or to better model results) + # the 18m figure is as a result of minimising pixel-wise errors between this model + # and that used in Done et al. 2020 with a physical boundary layer + GRADIENT_LEVEL_METRES = 18 + SURFACE_LEVEL_METRES = 10 + + logging.info("Calculating downscaling factors from surface roughness raster") + + # surface roughness raster for downscaling winds with + try: + surface_roughness_raster: xr.DataArray = rioxarray.open_rasterio(snakemake.input.surface_roughness) + except rasterio.errors.RasterioIOError: + logging.info("Surface roughness raster is empty, writing empty downscaling factors matrix.") + + # assure (empty) files exist + np.save(snakemake.output.downscale_factors, np.array([])) + os.system(f"touch {snakemake.output.downscale_factors_plot}") + + sys.exit(0) + + # (1, y, x) where 1 is the number of surface roughness bands + # select the only value in the band dimension + surface_roughness: np.ndarray = surface_roughness_raster[0].values + + # calculate factors to scale wind speeds from gradient-level to surface level, + # taking into account surface roughness as defined by the raster + downscaling_factors = power_law_scale_factors( + surface_roughness, + SURFACE_LEVEL_METRES, + GRADIENT_LEVEL_METRES + ) + + logging.info("Saving downscaling factors to disk") + np.save(snakemake.output.downscale_factors, downscaling_factors) + + logging.info("Mapping downscaling factors and saving to disk") + plot_downscale_factors( + downscaling_factors, + "Wind downscaling factors", + snakemake.output.downscale_factors_plot + ) + + logging.info("Done estimating downscaling factors") \ No newline at end of file diff --git a/workflow/scripts/preprocess/parse_IBTrACS.py b/workflow/scripts/preprocess/parse_IBTrACS.py index 38696513..caf94f54 100644 --- a/workflow/scripts/preprocess/parse_IBTrACS.py +++ b/workflow/scripts/preprocess/parse_IBTrACS.py @@ -11,7 +11,7 @@ import numpy as np import pandas as pd -from open_gira.io import STORM_CSV_SCHEMA +from parse_STORM import STORM_CSV_SCHEMA """ diff --git a/workflow/scripts/preprocess/parse_IRIS.py b/workflow/scripts/preprocess/parse_IRIS.py index 65969620..17911719 100644 --- a/workflow/scripts/preprocess/parse_IRIS.py +++ b/workflow/scripts/preprocess/parse_IRIS.py @@ -40,9 +40,10 @@ csv_dir = snakemake.input.csv_dir parquet_path = snakemake.output.parquet + sample = snakemake.wildcards.SAMPLE data = [] - for path in tqdm(natural_sort(glob(f"{csv_dir}/*.txt"))): + for path in tqdm(natural_sort(glob(f"{csv_dir}/*_1000Y_n{sample}.txt"))): df = pd.read_csv(path, header=1, sep=r"\s+", names=IRIS_CSV_SCHEMA.keys(), dtype=IRIS_CSV_SCHEMA) df = df.rename( diff --git a/workflow/scripts/preprocess/parse_STORM.py b/workflow/scripts/preprocess/parse_STORM.py index daeb93b6..3b14d741 100644 --- a/workflow/scripts/preprocess/parse_STORM.py +++ b/workflow/scripts/preprocess/parse_STORM.py @@ -48,17 +48,14 @@ csv_dir = snakemake.input.csv_dir parquet_path = snakemake.output.parquet + sample = snakemake.wildcards.SAMPLE data = [] - for path in tqdm(natural_sort(glob(f"{csv_dir}/*.csv"))): + # loop over basins for given sample, processing tracks into a common format + for path in tqdm(natural_sort(glob(f"{csv_dir}/*_1000_YEARS_{sample}*.csv"))): df = pd.read_csv(path, names=STORM_CSV_SCHEMA.keys(), dtype=STORM_CSV_SCHEMA) - # example paths containing sample number: - # STORM_DATA_HadGEM3-GC31-HM_WP_1000_YEARS_9_IBTRACSDELTA.csv - # STORM_DATA_IBTRACS_EP_1000_YEARS_0.csv - sample, = re.search(r"1000_YEARS_([\d])", os.path.basename(path)).groups() - df["sample"] = int(sample) # this gives us 10,000 years of data rather than 10 * 10,000 years