Fix bugs in make_state_vector.py

[laestrada]

Name and Institution (Required)

Name: Lucas Estrada Institution: Harvard ACMG

Describe the update

This PR addresses the following bugs in make_state_vector.py:

• an offset was mistakenly being applied in non 0.25 degree resolution grids
• faulty if statement caused land threshold to always evaluated to false
• In some cases labels would not be applied to buffer cells, causing gaps in the statevector labels resulting in perturbation simulations that would not be perturbed at all. This seems to only occur when using custom state vectors from a shapefile

To fix the 3rd bug I updated the buffer cell creation code and aligned them in src/utilities/make_state_vector_file.py and src/notebooks/statevector_from_shapefile.ipynb

@eastjames, tagging you and issue #181

[eastjames]

Thanks, Lucas. Just a couple questions on the first 2 items.

For bullet 1 it looks good. @nicholasbalasus said in #181 that updating the HEMCO standalone outputs would get rid of the offset issue completely, so this may not be needed after that happens.

For bullet 2, this assumes that land_threshold is always set to a numerical value. Is that the expected behavior?

[laestrada]

Hi @eastjamesm,

1. Yes, once the SA output is updated we can remove that line completely.
2. Yes, land_threshold is set by LandThreshold in the config file and will always be a value between 0 and 1

[eastjames]

Thanks, good to go.

[nicholasbalasus]

Hey Lucas, if I remember correctly the 0.5° x 0.625° output is offset by 0.0625°.

[nicholasbalasus]

This is untested and might not work with floating point issues, but maybe:

if np.abs(lc.lon.values - hd.lon.values).max() == 0.03125:
         hd["lon"] = hd["lon"] - 0.03125
elif np.abs(lc.lon.values - hd.lon.values).max() == 0.0625:
         hd["lon"] = hd["lon"] - 0.0625

[eastjames]

Thanks Nick, you're right, assuming the difference is 0.03125 isn't flexible for other grids. What about this for flexibility to avoid enumerating each case, or are those the only two cases? Also untested for floating point issues.

diff_threshold = 0.001
lon_diff = np.abs(lc.lon.values - hd.lon.values).max()
if lon_diff > diff_threshold:
    hd["lon"] = hd["lon"] - lon_diff

[nicholasbalasus]

Looks good to me!

[laestrada]

@eastjames and @nicholasbalasus,

I just did a quick check if the above would work and it looks like the IMI uses the nested-region specific landcover file, but the global file for emissions, so np.abs(lc.lon.values - hd.lon.values).max() yields an error due to differences in shapes. We could update it based on the bounds of the domain, but given that we are planning to update the SA output files soon anyways I am going to update the lines to:

if config["Res"] == "0.25x0.3125":
         hd["lon"] = hd["lon"] - 0.03125
elif config["Res"] == "0.5x0.625":
         hd["lon"] = hd["lon"] - 0.0625

and add a TODO reminding us to remove it when the diagnostics file is next updated

Updated

[eastjames]

Thanks Lucas. I tested with global so I didn’t catch the error with different array shapes