openghg / openghg_inversions

University of Bristol Atmospheric Chemistry Research Group RHIME Inversion code (with openghg dependency)
MIT License
5 stars 0 forks source link

postprocessing not working for 2-year inversions #223

Open hdelongueville opened 2 weeks ago

hdelongueville commented 2 weeks ago

When running inversions with a 2-year inversion period, RHIME encounters issues during post-processing. It would be nice to have the flexibility to run inversions with any inversion period.

Traceback (most recent call last): File "/user/home/qq24644/openghg_inversions/openghg_inversions/hbmcmc/run_hbmcmc.py", line 233, in mcmc_function(param) File "/user/home/qq24644/openghg_inversions/openghg_inversions/hbmcmc/hbmcmc.py", line 550, in fixedbasisMCMC out = mcmc.inferpymc_postprocessouts(post_process_args) ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ File "/user/home/qq24644/openghg_inversions/openghg_inversions/hbmcmc/inversion_pymc.py", line 721, in inferpymc_postprocessouts apriori_flux += flux_array_all[:, :, m] * np.sum(allmonths == m) / len(allmonths)


IndexError: index 2 is out of bounds for axis 2 with size 2
hdelongueville commented 1 week ago

It looks like the issue comes from inferpymc_postprocessouts, where if it assumes flux prior is monthly if not annual.

    if flux_array_all.shape[2] == 1:
        print("\nAssuming flux prior is annual and extracting first index of flux array.")
        apriori_flux = flux_array_all[:, :, 0]
    else:
        print("\nAssuming flux prior is monthly.")
        print(f"Extracting weighted average flux prior from {start_date} to {end_date}")
        allmonths = pd.date_range(start_date, end_date).month[:-1].values
        allmonths -= 1  # to align with zero indexed array

        apriori_flux = np.zeros_like(flux_array_all[:, :, 0])

        # calculate the weighted average flux across the whole inversion period
        for m in np.unique(allmonths):
            apriori_flux += flux_array_all[:, :, m] * np.sum(allmonths == m) / len(allmonths)
hdelongueville commented 1 week ago

@brendan-m-murphy I know you are working on updating the post-processing part. Do you have an opinion on the best way to solve this without interrupting your work (or maybe this will not impact what you are doing)?