From e236518dd5eca73a029ee345839d83d915771374 Mon Sep 17 00:00:00 2001 From: Lucas A Estrada <63303345+laestrada@users.noreply.github.com> Date: Wed, 4 Oct 2023 10:04:01 -0400 Subject: [PATCH] Feature/optimize domain edges (#169) This update implements boundary condition optimization as part of the inversion, making use of the new boundary condition perturbation options from @nicholasbalasus in geoschem 14.2.1. This is done by effectively increasing the number of state vector elements by 4, and performing 4 additional perturbation simulations. Notably, we do the inversion calculations for these final 4 elements in concentration space instead of calculating a correction factor. The calculated increase/decrease in ppb is then applied in the posterior simulation. These changes have been tested using our standard 1-week test case, but the differences are minimal when using buffer cells and BC optimization. Additional testing was done for 1-month over Tasmania to see the effect of having no buffer cells for the domain with and without BC optimization. --- config.yml | 7 ++ .../getting-started/imi-config-file.rst | 6 ++ envs/Harvard-Cannon/config.harvard-cannon.yml | 7 ++ resources/containers/container_config.yml | 7 ++ src/components/jacobian_component/jacobian.sh | 39 +++++++- src/components/kalman_component/kalman.sh | 2 +- .../posterior_component/posterior.sh | 26 +++++- src/components/setup_component/setup.sh | 2 +- .../statevector_component/statevector.sh | 2 +- src/inversion_scripts/calc_sensi.py | 36 +++++++- src/inversion_scripts/invert.py | 26 +++++- src/inversion_scripts/run_inversion.sh | 19 ++-- src/notebooks/visualization_notebook.ipynb | 88 ++++++------------- src/utilities/common.sh | 10 ++- src/utilities/sanitize_input_yaml.py | 3 +- 15 files changed, 194 insertions(+), 86 deletions(-) diff --git a/config.yml b/config.yml index 32d993c8..31a38c36 100644 --- a/config.yml +++ b/config.yml @@ -53,6 +53,7 @@ nBufferClusters: 8 BufferDeg: 5 LandThreshold: 0.25 OffshoreEmisThreshold: 0 +OptimizeBCs: false ## Point source datasets ## Used for visualizations and state vector clustering @@ -71,7 +72,10 @@ StateVectorFile: "/home/ubuntu/integrated_methane_inversion/resources/statevecto ShapeFile: "/home/ubuntu/integrated_methane_inversion/resources/shapefiles/PermianBasin_Extent_201712.shp" ## Inversion +## Note PriorError is a relative percentage +## and PriorErrorBCs is in ppb PriorError: 0.5 +PriorErrorBCs: 10.0 ObsError: 15 Gamma: 1.0 PrecomputedJacobian: false @@ -125,7 +129,10 @@ SchedulerPartition: "debug" ##-------------------------------------------------------------------- ## Jacobian settings +## Note PerturbValue is a relative scalefactor and +## PerturbValueBCs is in ppb PerturbValue: 1.5 +PerturbValueBCs: 10.0 ## Apply scale factors from a previous inversion? UseEmisSF: false diff --git a/docs/source/getting-started/imi-config-file.rst b/docs/source/getting-started/imi-config-file.rst index 2e8f974a..e9e27d08 100644 --- a/docs/source/getting-started/imi-config-file.rst +++ b/docs/source/getting-started/imi-config-file.rst @@ -85,6 +85,8 @@ State vector - Land-cover fraction below which to exclude GEOS-Chem grid cells from the state vector when creating the state vector file. Default value is ``0.25``. * - ``OffshoreEmisThreshold`` - Offshore GEOS-Chem grid cells with oil/gas emissions above this threshold will be included in the state vector. Default value is ``0``. + * - ``OptimizeBCs`` + - Boolean to optimize boundary conditions during the inversion. Must also include ``PerturbValueBCs`` and ``PriorErrorBCs`` Default value is ``false``. Clustering Options ^^^^^^^^^^^^^^^^^^ @@ -130,6 +132,8 @@ Inversion * - ``PriorError`` - Error in the prior estimates (1-sigma; relative). Default is ``0.5`` (50%) error. + * - ``PriorErrorBCs`` + - Error in the prior estimates (using ppb). Default is ``10`` ppb error. * - ``ObsError`` - Observational error (1-sigma; absolute; ppb). Default value is ``15`` ppb error. * - ``Gamma`` @@ -232,6 +236,8 @@ These settings are intended for advanced users who wish to modify additional GEO * - ``PerturbValue`` - Value to perturb emissions by in each sensitivity simulation. Default value is ``1.5``. + * - ``PerturbValueBCs`` + - Number of ppb to perturb emissions by for domain edges (North, South, East, West) if using `OptimizeBCs`. Default value is ``10.0`` ppb. * - ``UseEmisSF`` - Boolean to apply emissions scale factors derived from a previous inversion. This file should be provided as a netCDF file and specified in HEMCO_Config.rc. Default value is ``false``. * - ``UseOHSF`` diff --git a/envs/Harvard-Cannon/config.harvard-cannon.yml b/envs/Harvard-Cannon/config.harvard-cannon.yml index a612c21b..7bb79e41 100644 --- a/envs/Harvard-Cannon/config.harvard-cannon.yml +++ b/envs/Harvard-Cannon/config.harvard-cannon.yml @@ -53,6 +53,7 @@ nBufferClusters: 8 BufferDeg: 5 LandThreshold: 0.25 OffshoreEmisThreshold: 0 +OptimizeBCs: false ## Point source datasets ## Used for visualizations and state vector clustering @@ -71,7 +72,10 @@ StateVectorFile: "/n/holyscratch01/jacob_lab/msulprizio/Test_IMI/StateVector.nc" ShapeFile: "resources/shapefiles/PermianBasin_Extent_201712.shp" ## Inversion +## Note PriorError is a relative percentage +## and PriorErrorBCs is in ppb PriorError: 0.5 +PriorErrorBCs: 10.0 ObsError: 15 Gamma: 1.0 PrecomputedJacobian: false @@ -114,7 +118,10 @@ SchedulerPartition: "seas_compute,huce_cascade,huce_intel" ##-------------------------------------------------------------------- ## Jacobian settings +## Note PerturbValue is a relative scalefactor and +## PerturbValueBCs is in ppb PerturbValue: 1.5 +PerturbValueBCs: 10.0 ## Apply scale factors from a previous inversion? UseEmisSF: false diff --git a/resources/containers/container_config.yml b/resources/containers/container_config.yml index 6bf5be2d..24ed0815 100644 --- a/resources/containers/container_config.yml +++ b/resources/containers/container_config.yml @@ -53,6 +53,7 @@ nBufferClusters: 8 BufferDeg: 5 LandThreshold: 0.25 OffshoreEmisThreshold: 0 +OptimizeBCs: false ## Point source datasets ## Used for visualizations and state vector clustering @@ -71,7 +72,10 @@ StateVectorFile: "/home/al2/integrated_methane_inversion/resources/statevectors/ ShapeFile: "/home/al2/integrated_methane_inversion/resources/shapefiles/PermianBasin_Extent_201712.shp" ## Inversion +## Note PriorError is a relative percentage +## and PriorErrorBCs is in ppb PriorError: 0.5 +PriorErrorBCs: 10.0 ObsError: 15 Gamma: 1.0 PrecomputedJacobian: false @@ -125,7 +129,10 @@ SchedulerPartition: "debug" ##-------------------------------------------------------------------- ## Jacobian settings +## Note PerturbValue is a relative scalefactor and +## PerturbValueBCs is in ppb PerturbValue: 1.5 +PerturbValueBCs: 10.0 ## Apply scale factors from a previous inversion? UseEmisSF: false diff --git a/src/components/jacobian_component/jacobian.sh b/src/components/jacobian_component/jacobian.sh index 9eb5339e..a08ae2b4 100644 --- a/src/components/jacobian_component/jacobian.sh +++ b/src/components/jacobian_component/jacobian.sh @@ -3,6 +3,7 @@ # Functions available in this file include: # - setup_jacobian # - run_jacobian +# - generate_BC_perturb_values # Description: Setup jacobian run directory # Usage: @@ -81,9 +82,26 @@ setup_jacobian() { fi fi - # Update settings in geoschem_config.yml - sed -i -e "s|emission_perturbation: 1.0|emission_perturbation: ${PerturbValue}|g" \ - -e "s|state_vector_element_number: 0|state_vector_element_number: ${xUSE}|g" geoschem_config.yml + # Update settings in geoschem_config.yml except for the base run + if [ $x -ne 0 ]; then + sed -i -e "s|emission_perturbation: 1.0|emission_perturbation: ${PerturbValue}|g" \ + -e "s|state_vector_element_number: 0|state_vector_element_number: ${xUSE}|g" geoschem_config.yml + fi + + # BC optimization setup + if "$OptimizeBCs"; then + bcThreshold=$(($nElements - 4)) + # The last four state vector elements are reserved for BC optimization of NSEW + # domain edges. If the current state vector element is one of these, then + # turn on BC optimization for the corresponding edge and revert emission perturbation + if [ $x -gt $bcThreshold ]; then + PerturbBCValues=$(generate_BC_perturb_values $bcThreshold $x $PerturbValueBCs) + sed -i -e "s|CH4_boundary_condition_ppb_increase_NSEW:.*|CH4_boundary_condition_ppb_increase_NSEW: ${PerturbBCValues}|g" \ + -e "s|perturb_CH4_boundary_conditions: false|perturb_CH4_boundary_conditions: true|g" \ + -e "s|emission_perturbation: ${PerturbValue}|emission_perturbation: 1.0|g" \ + -e "s|state_vector_element_number: ${xUSE}|state_vector_element_number: 0|g" geoschem_config.yml + fi + fi # Update settings in HISTORY.rc # Only save out hourly pressure fields to daily files for base run @@ -185,3 +203,18 @@ run_jacobian() { printf "Got Jacobian scale factors\n" fi } + +# Description: Print perturbation string for BC optimization +# based on the current state vector element +# Returns [float, float, float, float] +# Usage: +# generate_BC_perturb_values +generate_BC_perturb_values() { + python -c "import sys;\ + bc_perturb = [0.0, 0.0, 0.0, 0.0];\ + bcThreshold = int(sys.argv[1]) + 1;\ + element = int(sys.argv[2]);\ + pert_index = element % bcThreshold;\ + bc_perturb[pert_index] = float(sys.argv[3]);\ + print(bc_perturb)" $1 $2 $3 +} diff --git a/src/components/kalman_component/kalman.sh b/src/components/kalman_component/kalman.sh index 5f774b4a..f40562cc 100644 --- a/src/components/kalman_component/kalman.sh +++ b/src/components/kalman_component/kalman.sh @@ -30,7 +30,7 @@ setup_kf() { mkdir -p ${RunDirs}/archive_sf # Number of state vector elements - nElements=$(ncmax StateVector ${StateVectorFile}) + nElements=$(ncmax StateVector ${StateVectorFile} ${OptimizeBCs}) } # Description: Run Kalman filter inversions diff --git a/src/components/posterior_component/posterior.sh b/src/components/posterior_component/posterior.sh index c850c9cd..c6de5e02 100644 --- a/src/components/posterior_component/posterior.sh +++ b/src/components/posterior_component/posterior.sh @@ -100,6 +100,21 @@ run_posterior() { source ${GEOSChemEnv} fi + if "$OptimizeBCs"; then + if "$KalmanMode"; then + inv_result_path="${RunDirs}/kf_inversions/period${period_i}/inversion_result.nc" + else + inv_result_path="${RunDirs}/inversion/inversion_result.nc" + fi + # set BC optimal delta values + PerturbBCValues=$(generate_optimized_BC_values $inv_result_path) + # add BC optimization delta to boundary condition edges + sed -i -e "s|CH4_boundary_condition_ppb_increase_NSEW:.*|CH4_boundary_condition_ppb_increase_NSEW: ${PerturbBCValues}|g" \ + -e "s|perturb_CH4_boundary_conditions: false|perturb_CH4_boundary_conditions: true|g" geoschem_config.yml + + printf "\n=== BC OPTIMIZATION: BC optimized perturbation values for NSEW set to: ${PerturbBCValues} ===\n" + fi + # Submit job to job scheduler printf "\n=== SUBMITTING POSTERIOR SIMULATION ===\n" sbatch --mem $SimulationMemory \ @@ -147,7 +162,7 @@ run_posterior() { LonMaxInvDomain=$(ncmax lon ${RunDirs}/StateVector.nc) LatMinInvDomain=$(ncmin lat ${RunDirs}/StateVector.nc) LatMaxInvDomain=$(ncmax lat ${RunDirs}/StateVector.nc) - nElements=$(ncmax StateVector ${RunDirs}/StateVector.nc) + nElements=$(ncmax StateVector ${RunDirs}/StateVector.nc ${OptimizeBCs}) FetchTROPOMI="False" isPost="True" buildJacobian="False" @@ -160,3 +175,12 @@ run_posterior() { # convert vizualization notebooks to html run_notebooks } + +# Description: Generates the updated NSEW perturbation to apply to domain edge BCs +# Usage: +# generate_optimized_BC_values +generate_optimized_BC_values() { + python -c "import sys; import xarray;\ + xhat = xarray.open_dataset(sys.argv[1])['xhat'].values[-4:];\ + print(xhat.tolist())" $1 +} \ No newline at end of file diff --git a/src/components/setup_component/setup.sh b/src/components/setup_component/setup.sh index 7d7e7dbe..45d13461 100644 --- a/src/components/setup_component/setup.sh +++ b/src/components/setup_component/setup.sh @@ -161,7 +161,7 @@ setup_imi() { fi # Determine number of elements in state vector file - nElements=$(ncmax StateVector ${RunDirs}/StateVector.nc) + nElements=$(ncmax StateVector ${RunDirs}/StateVector.nc ${OptimizeBCs}) printf "\nNumber of state vector elements in this inversion = ${nElements}\n\n" # Define inversion domain lat/lon bounds diff --git a/src/components/statevector_component/statevector.sh b/src/components/statevector_component/statevector.sh index c0522697..1912e093 100644 --- a/src/components/statevector_component/statevector.sh +++ b/src/components/statevector_component/statevector.sh @@ -99,7 +99,7 @@ reduce_dimension() { mkdir -p ${RunDirs}/archive_sv cp $state_vector_path ${RunDirs}/archive_sv/StateVector_${period_i}.nc fi - nElements=$(ncmax StateVector ${RunDirs}/StateVector.nc) + nElements=$(ncmax StateVector ${RunDirs}/StateVector.nc ${OptimizeBCs}) printf "\nNumber of state vector elements in this inversion = ${nElements}\n\n" printf "\n=== DONE REDUCING DIMENSION OF STATE VECTOR FILE ===\n" } \ No newline at end of file diff --git a/src/inversion_scripts/calc_sensi.py b/src/inversion_scripts/calc_sensi.py index 27d78bac..e5530160 100644 --- a/src/inversion_scripts/calc_sensi.py +++ b/src/inversion_scripts/calc_sensi.py @@ -15,9 +15,29 @@ def zero_pad_num(n): nstr = "0" + nstr return nstr +def test_GC_output_for_BC_perturbations(e, nelements, sensitivities): + + """ + Ensures that CH4 boundary condition perturbation in GEOS-Chem is working as intended + sensitivities = (pert-base)/perturbationBC which should equal 1e-9 inside the perturbed borders + example: the north boundary is perturbed by 10 ppb + pert-base=10e-9 mol/mol in the 3 grid cells that have been perturbed + perturbationBC=10 ppb + sensitivities = (pert-base)/perturbationBC = 1e-9 + """ + + if e == (nelements - 4): # North boundary + check = np.mean(sensitivities[:,-3:,3:-3]) + elif e == (nelements - 3): # South boundary + check = np.mean(sensitivities[:,0:3,3:-3]) + elif e == (nelements - 2): # East boundary + check = np.mean(sensitivities[:,:,-3:]) + elif e == (nelements - 1): # West boundary + check = np.mean(sensitivities[:,:,0:3]) + assert abs(check - 1e-9) < 1e-11, f"GC CH4 perturb not working... perturbation is off by {abs(check - 1e-9)} mol/mol/ppb" def calc_sensi( - nelements, perturbation, startday, endday, run_dirs_pth, run_name, sensi_save_pth + nelements, perturbation, startday, endday, run_dirs_pth, run_name, sensi_save_pth, perturbationBC ): """ Loops over output data from GEOS-Chem perturbation simulations to compute sensitivities @@ -25,12 +45,13 @@ def calc_sensi( Arguments nelements [int] : Number of state vector elements - perturbation [float] : Size of perturbation (e.g., 0.5) + perturbation [float] : Size of perturbation (e.g., 1.5) startday [str] : First day of inversion period; formatted YYYYMMDD endday [str] : Last day of inversion period; formatted YYYYMMDD run_dirs_pth [str] : Path to directory containing GC Jacobian run directories run_name [str] : Simulation run name; e.g. 'CH4_Jacobian' sensi_save_pth [str] : Path to save the sensitivity data + perturbationBC [float] : Size of BC perturbation in ppb (eg. 10.0) Resulting 'sensi' files look like: @@ -62,6 +83,8 @@ def calc_sensi( sensi[element,:,:,:] = sens save sensi as netcdf with appropriate coordinate variables """ + # subtract by 1 because here we assume .5 is a +50% perturbation + perturbation = perturbation - 1 # Make date range days = [] @@ -106,7 +129,12 @@ def process(h): # Get the data for the current hour pert = pert_data["SpeciesConcVV_CH4"][h, :, :, :] # Compute and store the sensitivities - sensitivities = (pert.values - base.values) / perturbation + if (perturbationBC is not None) and (e >= (nelements-4)): + sensitivities = (pert.values - base.values) / perturbationBC + if h != 0: # because we take the first hour on the first day from spinup + test_GC_output_for_BC_perturbations(e, nelements, sensitivities) + else: + sensitivities = (pert.values - base.values) / perturbation sensi[e, :, :, :] = sensitivities # Save sensi as netcdf with appropriate coordinate variables sensi = xr.DataArray( @@ -140,6 +168,7 @@ def process(h): run_dirs_pth = sys.argv[5] run_name = sys.argv[6] sensi_save_pth = sys.argv[7] + perturbationBC = float(sys.argv[8]) if len(sys.argv) > 8 else None calc_sensi( nelements, @@ -149,4 +178,5 @@ def process(h): run_dirs_pth, run_name, sensi_save_pth, + perturbationBC, ) diff --git a/src/inversion_scripts/invert.py b/src/inversion_scripts/invert.py index 7652486a..dd3302f0 100644 --- a/src/inversion_scripts/invert.py +++ b/src/inversion_scripts/invert.py @@ -18,6 +18,7 @@ def do_inversion( gamma=0.25, res="0.25x0.3125", jacobian_sf=None, + prior_err_bc=None, ): """ After running jacobian.py, use this script to perform the inversion and save out results. @@ -200,17 +201,35 @@ def do_inversion( # Inverse of prior error covariance matrix, inv(S_a) Sa_diag = np.zeros(n_elements) Sa_diag.fill(prior_err**2) + + # if optimizing boundary conditions, adjust for it in the inversion + bc_idx = n_elements + if prior_err_bc is not None: + # add prior error for BCs + # as the last 4 elements of the diagonal + Sa_diag[-4:] = prior_err_bc**2 + bc_idx -= 4 + inv_Sa = np.diag(1 / Sa_diag) # Inverse of prior error covariance matrix # Solve for posterior scale factors xhat ratio = np.linalg.inv(gamma * KTinvSoK + inv_Sa) @ (gamma * KTinvSoyKxA) - xhat = 1 + ratio + + # update scale factors by 1 to match what geoschem expects + # Note: if optimizing BCs, the last 4 elements are in concentration + # space, so we do not need to add 1 + # xhat = 1 + ratio + xhat = ratio.copy() + xhat[:bc_idx] += 1 # Posterior error covariance matrix S_post = np.linalg.inv(gamma * KTinvSoK + inv_Sa) # Averaging kernel matrix A = np.identity(n_elements) - S_post @ inv_Sa + + # Print some statistics + print("Min:", xhat[:bc_idx].min(), "Mean:", xhat[:bc_idx].mean(), "Max", xhat[:bc_idx].max()) return xhat, ratio, KTinvSoK, KTinvSoyKxA, S_post, A @@ -253,6 +272,7 @@ def calculate_superobservation_error(sO, p): gamma = float(sys.argv[10]) res = sys.argv[11] jacobian_sf = sys.argv[12] + prior_err_BC = float(sys.argv[13]) if len(sys.argv) > 13 else None # Reformat Jacobian scale factor input if jacobian_sf == "None": @@ -271,6 +291,7 @@ def calculate_superobservation_error(sO, p): gamma, res, jacobian_sf, + prior_err_BC, ) xhat = out[0] ratio = out[1] @@ -279,9 +300,6 @@ def calculate_superobservation_error(sO, p): S_post = out[4] A = out[5] - # Print some statistics - print("Min:", xhat.min(), "Mean:", xhat.mean(), "Max", xhat.max()) - # Save results dataset = Dataset(output_path, "w", format="NETCDF4_CLASSIC") nvar = dataset.createDimension("nvar", n_elements) diff --git a/src/inversion_scripts/run_inversion.sh b/src/inversion_scripts/run_inversion.sh index cf7ee4ad..1bbf439b 100755 --- a/src/inversion_scripts/run_inversion.sh +++ b/src/inversion_scripts/run_inversion.sh @@ -103,14 +103,14 @@ printf "DONE -- postproc_diags.py\n\n" #======================================================================= if ! "$PrecomputedJacobian"; then - - # 50% perturbation - Perturbation=0.5 - + python_args=(calc_sensi.py $nElements $PerturbValue $StartDate $EndDate $JacobianRunsDir $RunName $sensiCache) + # add an argument to calc_sensi.py if optimizing BCs + if "$OptimizeBCs"; then + python_args+=($PerturbValueBCs) + fi printf "Calling calc_sensi.py\n" - python calc_sensi.py $nElements $Perturbation $StartDate $EndDate $JacobianRunsDir $RunName $sensiCache; wait + python "${python_args[@]}"; wait printf "DONE -- calc_sensi.py\n\n" - fi #======================================================================= @@ -158,8 +158,13 @@ fi posteriorSF="./inversion_result.nc" +python_args=(invert.py $nElements $JacobianDir $posteriorSF $LonMinInvDomain $LonMaxInvDomain $LatMinInvDomain $LatMaxInvDomain $PriorError $ObsError $Gamma $Res $jacobian_sf) +# add an argument to calc_sensi.py if optimizing BCs +if "$OptimizeBCs"; then + python_args+=($PriorErrorBCs) +fi printf "Calling invert.py\n" -python invert.py $nElements $JacobianDir $posteriorSF $LonMinInvDomain $LonMaxInvDomain $LatMinInvDomain $LatMaxInvDomain $PriorError $ObsError $Gamma $Res $jacobian_sf; wait +python "${python_args[@]}"; wait printf "DONE -- invert.py\n\n" #======================================================================= diff --git a/src/notebooks/visualization_notebook.ipynb b/src/notebooks/visualization_notebook.ipynb index 4d1e7588..14b466ff 100644 --- a/src/notebooks/visualization_notebook.ipynb +++ b/src/notebooks/visualization_notebook.ipynb @@ -259,83 +259,49 @@ ")" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Sectoral Emissions" + ] + }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ - "# Calculate posterior sectoral emissions\n", - "prior = xr.load_dataset(prior_pth)\n", - "posterior = prior * scale\n", - "\n", - "emission_types = {\n", - " \"SoilPost emissions\": sum_total_emissions(\n", - " posterior[\"EmisCH4_SoilAbsorb\"].isel(time=0), areas, mask\n", - " ),\n", - " \"Termites emissions\": sum_total_emissions(\n", - " posterior[\"EmisCH4_Termites\"].isel(time=0), areas, mask\n", - " ),\n", - " \"Lakes emissions\": sum_total_emissions(\n", - " posterior[\"EmisCH4_Lakes\"].isel(time=0), areas, mask\n", - " ),\n", - " \"Seeps emissions\": sum_total_emissions(\n", - " posterior[\"EmisCH4_Seeps\"].isel(time=0), areas, mask\n", - " ),\n", - " \"Wetlands emissions\": sum_total_emissions(\n", - " posterior[\"EmisCH4_Wetlands\"].isel(time=0), areas, mask\n", - " ),\n", - " \"BiomassBurn emissions\": sum_total_emissions(\n", - " posterior[\"EmisCH4_BiomassBurn\"].isel(time=0), areas, mask\n", - " ),\n", - " \"OtherAnth emissions\": sum_total_emissions(\n", - " posterior[\"EmisCH4_OtherAnth\"].isel(time=0), areas, mask\n", - " ),\n", - " \"Rice emissions\": sum_total_emissions(\n", - " posterior[\"EmisCH4_Rice\"].isel(time=0), areas, mask\n", - " ),\n", - " \"Wastewater emissions\": sum_total_emissions(\n", - " posterior[\"EmisCH4_Wastewater\"].isel(time=0), areas, mask\n", - " ),\n", - " \"Landfills emissions\": sum_total_emissions(\n", - " posterior[\"EmisCH4_Landfills\"].isel(time=0), areas, mask\n", - " ),\n", - " \"Livestock emissions\": sum_total_emissions(\n", - " posterior[\"EmisCH4_Livestock\"].isel(time=0), areas, mask\n", - " ),\n", - " \"Coal emissions\": sum_total_emissions(\n", - " posterior[\"EmisCH4_Coal\"].isel(time=0), areas, mask\n", - " ),\n", - " \"Gas emissions\": sum_total_emissions(\n", - " posterior[\"EmisCH4_Gas\"].isel(time=0), areas, mask\n", - " ),\n", - "}\n", + "def my_autopct(\n", + " pct,\n", + "): # Only plots the data labels for sectoral emissions that are larger than 15%\n", + " return (\"%.1f%%\" % pct) if pct > 15 else \"\"\n", "\n", - "labels = []\n", - "values = []\n", + "# load in full sectors based on prior distribution\n", + "posterior_sectors = xr.load_dataset(prior_pth) * scale\n", "\n", - "for subemission in emission_types:\n", - " print(subemission, \"=\", emission_types[subemission], \"Tg/y\")\n", - " if emission_types[subemission] > 0:\n", - " labels.append(subemission)\n", - " values.append(emission_types[subemission])\n", + "# extract sector names\n", + "sector_list = [var for var in list(posterior_sectors.keys()) if \"EmisCH4\" in var]\n", + "sector_list.remove(\"EmisCH4_Total\")\n", "\n", + "# calculate total emissions for each sector and print\n", + "emission_types = {}\n", + "for sector in sector_list:\n", + " emission = sum_total_emissions(\n", + " posterior_sectors[sector].isel(time=0), areas, mask\n", + " )\n", + " print(f\"{sector} = {emission} Tg/yr\")\n", + " if emission > 0:\n", + " emission_types[sector] = emission\n", "\n", "title = plt.title(\"CH4 emissions by sector\")\n", "title.set_ha(\"center\")\n", "plt.gca().axis(\"equal\")\n", "\n", - "\n", - "def my_autopct(\n", - " pct,\n", - "): # Only plots the data labels for sectoral emissions that are larger than 15%\n", - " return (\"%.1f%%\" % pct) if pct > 15 else \"\"\n", - "\n", - "\n", - "pie = plt.pie(values, startangle=0, autopct=my_autopct)\n", + "pie = plt.pie(np.array(list(emission_types.values()))/total_posterior_emissions, startangle=0, autopct=my_autopct)\n", "plt.legend(\n", " pie[0],\n", - " labels,\n", + " emission_types.keys(),\n", " bbox_to_anchor=(1, 0.5),\n", " loc=\"center right\",\n", " fontsize=10,\n", diff --git a/src/utilities/common.sh b/src/utilities/common.sh index 4a9126ea..1e4dff4b 100644 --- a/src/utilities/common.sh +++ b/src/utilities/common.sh @@ -36,11 +36,15 @@ imi_failed() { # Description: Print max value of given variable in netCDF file # Returns int if only trailing zeros, float otherwise +# if (optional) 3rd argument is true, add 4 to the max value +# this is useful for adding in elements to account for +# optimization of BCs # Usage: -# ncmax +# ncmax ncmax() { - python -c "import sys; import xarray; \ - print('%g' % xarray.open_dataset(sys.argv[2])[sys.argv[1]].max())" $1 $2 + python -c "import sys; import xarray;\ + bc_offset = 4 if len(sys.argv) > 3 and sys.argv[3] == 'true' else 0;\ + print('%g' % (xarray.open_dataset(sys.argv[2])[sys.argv[1]].max()+bc_offset))" $1 $2 $3 } # Description: Print min value of given variable in netCDF file diff --git a/src/utilities/sanitize_input_yaml.py b/src/utilities/sanitize_input_yaml.py index f9e89baf..b050f471 100644 --- a/src/utilities/sanitize_input_yaml.py +++ b/src/utilities/sanitize_input_yaml.py @@ -110,6 +110,7 @@ "S3UploadPath", "S3UploadFiles", ] +conditional_dict["OptimizeBCs"] = ["PerturbValueBCs", "PriorErrorBCs"] def raise_error_message(var): """ @@ -123,7 +124,7 @@ def raise_error_message(var): + "https://imi.readthedocs.io/en/latest/getting-started/imi-config-file.html" ) raise ValueError(message) - + if __name__ == "__main__": config_path = sys.argv[1]