Compute bias instead of correlation in compare_salinity.py

[sloosvel]

Description

Update sea surface salinity diagnostic for CMUG deliverable D5.4

@portegam @kserradell I am opening this as a draft PR to have a first look at the updated results.

Here are the updated figures that are generated when computing the bias instead of the correlation, as indicated:

Please let me know if that is what is expected

• Closes Fixes in diagnostic script for recipe_sea_surface_salinity #2641
• Link to documentation:

---

Before you get started

• ☝ Create an issue to discuss what you are going to do

Checklist

It is the responsibility of the author to make sure the pull request is ready to review. The icons indicate whether the item will be subject to the 🛠 Technical or 🧪 Scientific review.

• 🛠 This pull request has a descriptive title
• 🛠 Code is written according to the code quality guidelines
• 🛠 Documentation is available
• 🛠 Tests run successfully
  - [ ] 🛠 The list of authors is up to date
  - [ ] 🛠 Any changed dependencies have been added or removed correctly
• 🛠 All checks below this pull request were successful

New or updated recipe/diagnostic

• 🧪 Recipe runs successfully
• 🧪 Recipe is well documented
• 🧪 Figure(s) and data look as expected from literature
• 🛠 Provenance information has been added

To help with the number of pull requests:

• 🙏 We kindly ask you to review two other open pull requests in this repository

[sloosvel]

You can double check the code changes in the Files changed tab

[portegam]

Thanks @sloosvel, this is very promising for a start. I see that some adjustments would be needed to get it in good shape.
These would be my first suggestions:

• Can you exclude the two mediterranean regions from the plot? They seem to be way off in terms of bias when compared with the other areas and the values they take are out of range in the radar plot.
• I've noticed that the scale of values in the radar plot is not showing. I think that it's because we also need to adapt plt.ylim. Can you change the line 122 to: plt.ylim(0, 10)? We'll probably need 1-2 further iterations to get it right.
• I might want to consider a different set of observations instead of the ESA-SEA-SURFACE-SALINITY_V1. After checking the code, I couldn't find the line where the observations are loaded. Is this something that is requested outside of the function compare_salinity.py?

[sloosvel]

Great, thanks for the feedback. I will get on the changes.

I might want to consider a different set of observations instead of the ESA-SEA-SURFACE-SALINITY_V1. After checking the code, I couldn't find the line where the observations are loaded. Is this something that is requested outside of the function compare_salinity.py?

Yes, the datasets are specified in the recipe under the dataset section:

https://github.com/ESMValGroup/ESMValTool/blob/main/esmvaltool/recipes/recipe_sea_surface_salinity.yml

ESMValTool/esmvaltool/recipes/recipe_sea_surface_salinity.yml

Line 47 in d8687d3

datasets:

[sloosvel]

Hi @portegam ,

Here you have the second batch of result, the plot limits were updated and the Mediterranean regions were removed:

[portegam]

Yes, the datasets are specified in the recipe under the dataset section:

https://github.com/ESMValGroup/ESMValTool/blob/main/esmvaltool/recipes/recipe_sea_surface_salinity.yml

Great, thanks. I've checked here and there is another supported dataset of sos for ESMValTool: that is WOA. Can you see if you can load it instead of the V1 from ESA?

[portegam]

Hi @portegam ,

Here you have the second batch of result, the plot limits were updated and the Mediterranean regions were removed:

Thanks @sloosvel. I think I understand now how the y-axis work.
I'd like to do now a test to see if it's possible to make a plot that looks good with biases that are not absolute (and can be either positive or negative). Could you redo the plots with?

line 101> bias_data = climat_ref.data - climat_data.data
line 119> plt.yticks([-0.5, 0.0,0.5], ["-0.5", "0", "0.5"],
line 122> plt.ylim(-1, 1)
line 149> caption = (f"Absolute bias comparison in different regions for "

If it looks good enough, it might be the final version of the script

[sloosvel]

Great, thanks. I've checked here and there is another supported dataset of sos for ESMValTool: that is WOA. Can you see if you can load it instead of the V1 from ESA?

Which version would you need? 2018 or 2013v2? Also it seems like it's a climatology value so I would only be able to plug it in the bias plots, but not the in timeseries ones.

[portegam]

Which version would you need? 2018 or 2013v2? Also it seems like it's a climatology value so I would only be able to plug it in the bias plots, but not the in timeseries ones.

If it's a climatology then I think that it's better not to use it, because it won't be computed for the same period that we have used (that is the overlap between the obs and the model). The only other product with salinity data (3D instead of only surface) that I've found is PHC. Can you check with versions are available and if they provide time varying data or just a climatology?

[sloosvel]

PHC also contains just one time point:

dimensions:
        time = UNLIMITED ; // (1 currently)
        lev = 33 ;
        lat = 180 ;
        lon = 360 ;
        bnds = 2 ;
...
:modeling_realm = "clim" ;
:comment = "The data are multiyear annual mean, but given Omon mip." ;

[portegam]

PHC also contains just one time point:

Thanks for checking @sloosvel.
I'm a bit surprised I couldn't find EN4, which is the most widely used dataset for ocean temperature and salinity, is not included in the list of https://docs.esmvaltool.org/en/latest/input.html.
Do you know if it's possible to load other datasets (like EN4) that are not in the list?

[sloosvel]

Do you know if it's possible to load other datasets (like EN4) that are not in the list?

I would need to add a cmoriser script first

[portegam]

I would need to add a cmoriser script first
It's good to know that it's possible. Ok, I'd suggest to skip it, and keep the two ESA products as it's currently done. Just focus on the other changes I suggested for the radar plot.

[sloosvel]

And the plots for the newest scale:

[portegam]

And the plots for the newest scale:
Great, I think this will work. Would it be difficult to add the recipe for a new type of radar plot. It would be very similar to the one used for the biases, but with a few line changes.

[sloosvel]

Sure, no problem! What type of plot do you further need?

[portegam]

It would be a similar plot, but instead of showing the bias (difference in mean states), it would show the ratio between the observed and the simulated standard deviation (calculated both over the overlap period). I will copy below the changes in the compare_salinity.py script that would need to be done to make that computation. The script for that could be called compare_salinity_std.py to differentiate it from the previous.

[portegam]

So the main changes to do would be:

• In lines 99 and 100, change the 'mean' operator called by climate_statistics to the standard deviation (I couldn't find how it's called in the package, but I guess that it's std)
• In line 101, instead of "bias_data = climat_ref.data - climat_data.data" it would be the equivalent to "ratio_std = climat_ref.data / climat_data.data" (note that now bias_data would need to be changed to ratio_std, and the bias variable should be changed too, e.g., to stdr)
  We would also need to adjust the y axis in the radar function, as the range of values is expected to be different. I would start with:
• Line 119: plt.yticks([0.25, 0.5, 1, 2, 4], ["0.25", "0.5", "1", "2","4"],
• Line 122: plt.ylim(0, 5)

[sloosvel]

Please find below the std ratio results:

[portegam]

Very nice plots, thanks @sloosvel. Do you know if there is an option to plot the y-axis on a logarithmic basis? It might improve the visibility of the results because the interpretation of the values that are greater/lower than 1 is non linear.

[sloosvel]

Yes, but the ticks change to logarithmic. Would that work for you?

[portegam]

Thanks @sloosvel. It might work. Can you test it for the range 10^-2 to 10^2? (with ticks on 10^-1, 10^0 and 10^1)?

[portegam]

I'd like also if you could repeat the plot in the linear (non-logaritmic) scale, using the same yticks "0, 0.25, 0.5, 1, 2, 4" and the range -2 to 5.

[sloosvel]

Here you have the results:

[portegam]

I don't have access to jasmin, I'm afraid. Do you mind uploading here just that figure?

[sloosvel]

Sure:

[portegam]

Thanks a lot @sloosvel

[sloosvel]

Hi @portegam, I did a small update to the documentation. Could you please check if that reads fine to you and if so just say that you give your scientific approval? Thanks in advance.

[sloosvel]

The technical approval was given in #1832 and the current updates do not affect the quality of the code, as all checks pass.

[portegam]

Hi @portegam, I did a small update to the documentation. Could you please check if that reads fine to you and if so just say that you give your scientific approval? Thanks in advance.
Sure, no problem. I'm into it.

[sloosvel]

Then you would need to go to the Files changed tab, click on the green Review changes button, select the Approve option and click on the green Submit review button, otherwise I won't be able to merge this. Thanks.

[portegam]

I added 3 comments with suggested changes in the text (nothing was wrong, I just added a bit more detail).

Then you would need to go to the Files changed tab, click on the green Review changes button, select the Approve option and click on the green Submit review button, otherwise I won't be able to merge this. Thanks.

I think I did it that way. Please let me know if it's not the case.

[portegam]

I also have one question. Javi's recipe only computed one type of radar plot, and with our changes we are computing two. Are both computed by the same script? Or do we have two separate scripts?

[sloosvel]

I also have one question. Javi's recipe only computed one type of radar plot, and with our changes we are computing two. Are both computed by the same script? Or do we have two separate scripts?

It's one script that loops over the mean and std_dev operators, in that way there is less code duplication.

[portegam]

It's one script that loops over the mean and std_dev operators, in that way there is less code duplication.

Excellent, it's much better than way. Thanks for implementing it that way, it's much cleaner

[portegam]

Hi @sloosvel, is there anything else needed from my side? I see that the required review appears to be still missing.

[sloosvel]

Actually I think I don't have permissions to merge so if anyone from @ESMValGroup/esmvaltool-coreteam could do it that would be great.

[zklaus]

At the moment, merging is blocked because there is no review.

[zklaus]

@portegam, are you a member of ESMValGroup/esmvaltool-developmentteam ?

[sloosvel]

@pepcos can you do us a favour and approve the changes?

[portegam]

I'm not, I'm just the BSC person that has been providing the scientific guidanceEl 6 may. 2022 15:19, Klaus Zimmermann ***@***.***> escribió: @portegam, are you a member of ESMValGroup/esmvaltool-developmentteam ? —Reply to this email directly, view it on GitHub, or unsubscribe.You are receiving this because you were mentioned.Message ID: ***@***.***>

[zklaus]

Ah, that's why Github didn't count your review. Was just wondering, thanks for clarifying.

[sloosvel]

Thanks!!