Zonal mean plots show striping when constant values everwhere

[lizziel]

Your name

Lizzie Lundgren

Your affiliation

Harvard University

What happened? What did you expect to happen?

I have noticed that using the latest dev results in strange zonal mean plots when plotting an array that is constant values everywhere. I am seeing this when I trouble-shoot a bug in GCHP passive tracer monthly mean diagnostics. The bug is that values are constant. The level plots look as expected, but the zonal mean looks like this:

What are the steps to reproduce the bug?

Make a zonal mean plot of GCHP passive tracer array (SpeciesConc collection variable SpeciesConcVV_PassiveTracer).

Please attach any relevant configuration and log files.

Steps to reproduce on Cannon: file1=~elundgren/testruns/gchp_14.4.0_test/gchp_merra2_TransportTracers/OutputDir/GEOSChem.SpeciesConc.20190101_0000z.nc4

Using ~elundgren/bin/compare_var.py:
compare_var.py $file1 'Monthly mean' $file1 'Monthly mean' SpeciesConcVV_PassiveTracer zonal 0 passive_tracer_zonal_monthly_mean

This will create file passive_tracer_zonal_monthly_mean.pdf which shows the issue using one dataset compared with itself.

What GCPy version were you using?

dev branch as of June 5, 2024 (a0172fc, Merge PR #327)

What environment were you running GCPy on?

Local cluster

What python environment were you using?

GCPy v1.4.1 environment

Will you be addressing this bug yourself?

Yes

Additional information

No response

[yantosca]

@lizziel: I ran this plot on your sample file using newer Python packages and got this output:

Which looks better except for the symmetric striping at about 40N and 40S latitude. I will try to see what is causing that.

[yantosca]

@lizziel: I think it may be related to the regridding to 1 x 1.25. I added this line into compare_zonal_mean:

        # ==============================================================
        # Assign data variables
        # ==============================================================
        ds_ref = ds_refs[ivar]
        ds_dev = ds_devs[ivar]
        ds_ref_cmp = ds_ref_cmps[ivar]
        ds_dev_cmp = ds_dev_cmps[ivar]
        frac_ds_ref_cmp = frac_ds_ref_cmps[ivar]
        frac_ds_dev_cmp = frac_ds_dev_cmps[ivar]

        print(f"{ivar}, {np.min(ds_ref_cmps[ivar])}, {np.max(ds_ref_cmps[ivar])}")

and re-ran the test case. The output was

0, 100.00063705444336, 100.00064468383789

so it could be some numerical drift in the regridding.

[yantosca]

@lizziel: I also put in this print statement

        # ==============================================================
        # Calculate zonal mean
        # ==============================================================
        # Ref
        if refgridtype == "ll":
            zm_ref = ds_ref.mean(dim="lon")
        else:
            zm_ref = ds_ref.mean(axis=2)

        # Dev
        if devgridtype == "ll":
            zm_dev = ds_dev.mean(dim="lon")
        else:
            zm_dev = ds_dev.mean(axis=2)
        # Comparison
        zm_dev_cmp = ds_dev_cmp.mean(axis=2)
        zm_ref_cmp = ds_ref_cmp.mean(axis=2)
        if diff_of_diffs:
            frac_zm_dev_cmp = frac_ds_dev_cmp.mean(axis=2)
            frac_zm_ref_cmp = frac_ds_ref_cmp.mean(axis=2)

        print(f"{np.min(zm_ref_cmp)}, {np.max(zm_ref_cmp)}")

which returns

100.00064076317682, 100.00064126060654

[yantosca]

I added a print statement to print out a latitudinal slice of the zonal mean data at the surface:

        # ==============================================================
        # Calculate zonal mean
        # ==============================================================
        # Ref
        print(type(ds_ref), type(ds_ref[0,0,0]))
        if refgridtype == "ll":
            zm_ref = ds_ref.mean(dim="lon")
        else:
            zm_ref = ds_ref.mean(axis=2)

        # Dev
        if devgridtype == "ll":
            zm_dev = ds_dev.mean(dim="lon")
        else:
            zm_dev = ds_dev.mean(axis=2)
        # Comparison
        zm_dev_cmp = ds_dev_cmp.mean(axis=2)
        zm_ref_cmp = ds_ref_cmp.mean(axis=2)
        if diff_of_diffs:
            frac_zm_dev_cmp = frac_ds_dev_cmp.mean(axis=2)
            frac_zm_ref_cmp = frac_ds_ref_cmp.mean(axis=2)

        print(zm_ref[0,:])

Which yields:

<class 'numpy.ndarray'> <class 'numpy.float64'>
[100.00064087 100.00064087 100.00064087 100.00064087 100.00064087
 100.00064087 100.00064087 100.00064087 100.00064087 100.00064087
 100.00064087 100.00064087 100.00064087 100.00064087 100.00064087
 100.00064087 100.00064087 100.00064087 100.00064087 100.00064087
 100.00064087 100.00064087 100.00064087 100.00064087 100.00064087
 100.00064087 100.00064087 100.00064087 100.00064087 100.00064087
 100.00064087 100.00064087 100.00064087 100.00064087 100.00064087
 100.00064087 100.00064087 100.00064087 100.00064087 100.00064087
 100.00064087 100.00064087 100.00064087 100.00064087 100.00064087
 100.00064126 100.00064102 100.00064098 100.00064088 100.00064089    <--- non-uniformity here
 100.00064095 100.00064107 100.00064076 100.00064082 100.00064086    <--- non-uniformity here
 100.00064092 100.00064087 100.00064087 100.00064087 100.00064087    <--- non-uniformity here
 100.00064087 100.00064087 100.00064087 100.00064087 100.00064087
 100.00064087 100.00064087 100.00064087 100.00064087 100.00064087
 100.00064087 100.00064087 100.00064087 100.00064087 100.00064087
 100.00064087 100.00064087 100.00064087 100.00064087 100.00064087
 100.00064087 100.00064087 100.00064087 100.00064087 100.00064087
 100.00064087 100.00064087 100.00064087 100.00064087 100.00064087
 100.00064087 100.00064087 100.00064087 100.00064087 100.00064087
 100.00064087 100.00064087 100.00064087 100.00064087 100.00064087
 100.00064087 100.00064087 100.00064087 100.00064087 100.00064087
 100.00064087 100.00064087 100.00064087 100.00064087 100.00064087
 100.00064087 100.00064087 100.00064087 100.00064087 100.00064087
 100.00064087 100.00064087 100.00064087 100.00064087 100.00064087
 100.00064087 100.00064087 100.00064087 100.00064087 100.00064087
 100.00064092 100.00064086 100.00064082 100.00064076 100.00064107  <--- non-uniformity here
 100.00064095 100.00064089 100.00064088 100.00064098 100.00064102  <--- non-uniformity here
 100.00064126 100.00064087 100.00064087 100.00064087 100.00064087  <--- non-uniformity here
 100.00064087 100.00064087 100.00064087 100.00064087 100.00064087
 100.00064087 100.00064087 100.00064087 100.00064087 100.00064087
 100.00064087 100.00064087 100.00064087 100.00064087 100.00064087
 100.00064087 100.00064087 100.00064087 100.00064087 100.00064087
 100.00064087 100.00064087 100.00064087 100.00064087 100.00064087
 100.00064087 100.00064087 100.00064087 100.00064087 100.00064087
 100.00064087 100.00064087 100.00064087 100.00064087 100.00064087
 100.00064087 100.00064087 100.00064087 100.00064087 100.00064087
 100.00064087]

Which is responsible for the striping that we see in the output.

Also tagging @sdeastham