Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix stats for climatology maps #977

Merged
merged 1 commit into from
Dec 11, 2023

Conversation

xylar
Copy link
Collaborator

@xylar xylar commented Dec 7, 2023

The area weights (and not just the fields themselves) need to be masked to exclude land.

Before this fix:
image
Note: mean of 12.36
But this time series of SST:
image

After the fix:
image

The area weights (and not just the fields themselves) need to
be masked to exclude land.
@xylar xylar self-assigned this Dec 7, 2023
@xylar xylar added the bug label Dec 7, 2023
@xylar
Copy link
Collaborator Author

xylar commented Dec 7, 2023

@golaz, I really though I had the stats fixed after #860 and #861 but clearly I messed up big time!

My only defense, and it's a weak one, is that it seems pretty odd that np.average() with a masked array doesn't apply the mask to the weights array. It seems pretty darn obvious that you wouldn't use a masked array if the mask values should just be treated as zeros for the purposes of the average. But I absolutely should have done more checks on the data before remapping.

@xylar
Copy link
Collaborator Author

xylar commented Dec 7, 2023

Testing

I ran the climatologyMapSST task with this fix:
https://web.lcrc.anl.gov/public/e3sm/diagnostic_output/ac.xylar/analysis/fix-climatology-map-stats/20231122.v3b02-Icos30.piControl.chrysalis/ts_0001-0300_climo_0251-0300/ocean/index.html#global_sst

I also computed the mean outside of MPAS-Analysis using both the MPAS data before remapping:

#!/usr/bin/env python

import xarray as xr

mesh_filename = '/lcrc/group/e3sm/public_html/inputdata/share/meshes/mpas/ocean/IcoswISC30E3r5.20231120.nc'
sst_filename = '/lcrc/group/e3sm2/ac.wlin/E3SMv3_dev/20231122.v3b02-Icos30.piControl.chrysalis/post/analysis/mpas_analysis/ts_0001-0300_climo_0251-0300/clim/mpas/avg/masked/sst_IcoswISC30E3r5/mpaso_ANN_025101_030012_climo.nc'
ds_mesh = xr.open_dataset(mesh_filename)
ds_sst = xr.open_dataset(sst_filename)

area = ds_mesh.areaCell

mean_sst = ((ds_sst.timeMonthly_avg_activeTracers_temperature * area).sum() /
            area.sum())

print(mean_sst.values)

The result is:
17.6265073782531

And the remapped data:

#!/usr/bin/env python

import xarray as xr


sst_filname = '/lcrc/group/e3sm2/ac.wlin/E3SMv3_dev/20231122.v3b02-Icos30.piControl.chrysalis/post/analysis/mpas_analysis/ts_0001-0300_climo_0251-0300/clim/mpas/avg/remapped/sst_IcoswISC30E3r5_to_0.5x0.5degree/mpaso_ANN_025101_030012_climo.nc'

ds_sst = xr.open_dataset(sst_filname)

area = ds_sst.area

sst = ds_sst.timeMonthly_avg_activeTracers_temperature

mask = sst.notnull()


mean_sst = ((sst * area).sum() / area.where(mask).sum())

print(mean_sst.values)

mean_sst = ((sst * area).sum() / area.sum())

print(mean_sst.values)

The results are:
17.67086178928455
when properly masking areas and
12.355902147205343
when areas are not masked.
So clearly the areas were not getting masked as they should have been.

Copy link
Collaborator

@milenaveneziani milenaveneziani left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Looks good to me. Thanks @xylar for the super quick fix!

@golaz
Copy link

golaz commented Dec 7, 2023

@xylar , thanks for fixing this so quickly. I agree that the behavior is confusing. For my own scripts, I typically use ma.average() when using masked arrays. I view ma as a layer on top of np. So ma functions will behave like their corresponding np versions if passed np arrays, but the reverse is not assured.

This bug also explains why SST RMSE metrics were always much smaller from MPAS-Analysis than from E3SM Diags. Something I had noticed a while ago, but did not dig into.

@xylar
Copy link
Collaborator Author

xylar commented Dec 7, 2023

@golaz, thanks for the advice on ma.average(). I didn't know that existed.

@milenaveneziani, I agree with your slack comment, I prefer the xarray options over np but this code is old enough that I must not have felt comfortable using them back then.

@xylar xylar removed the request for review from proteanplanet December 11, 2023 04:26
@xylar
Copy link
Collaborator Author

xylar commented Dec 11, 2023

Thanks @milenaveneziani and @golaz!

@proteanplanet, I think this has enough eyes on it already so I'll remove you as a reviewer. You have a busy week ahead of you!

@xylar xylar merged commit e7a1628 into MPAS-Dev:develop Dec 11, 2023
5 checks passed
@xylar xylar deleted the fix-climatology-map-stats branch December 11, 2023 04:28
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants