Add risk index outcome for navigation in ice-covered waters #968

torotorotaxi · 2023-07-11T08:15:31Z

This adds a new analysis task to compute the risk index outcome (rio) from sea-ice concentration and (floe) thickness, as outlined in the International Maritime Organization (IMO) document (https://www.imorules.com/GUID-2C1D86CB-5D58-490F-B4D4-46C057E1D102.html).

torotorotaxi · 2023-07-11T08:17:52Z

Plots need to be refined (mostly labels and titles). Probably some commits will need to be squashed before merging.

torotorotaxi · 2023-07-11T08:20:19Z

From Milena: @xylar, I helped Gennaro with this one. Hopefully, everything should be almost ready to go, but please let us know if you see anything amiss. Thanks. We'll post an example plot as soon as I can find it (we are in Italy right now, so same time zone as you).

torotorotaxi · 2023-07-11T08:25:33Z

Here is an example output for this task:
https://portal.nersc.gov/project/m1199/milena/E3SMv2.1B60to10rA02/Years96-100/sea_ice/index.html

@proteanplanet: it would also be great if you could take a look and comment here. It would be nice if we could get rid of the rio>30 noisy results (darkest blue), but the things the we tried didn't work.

milenaveneziani · 2023-07-12T08:02:08Z

@xylar, @proteanplanet, @stephenprice: in case you missed this since it came from a different username, I am pinging you again. No rush, just making sure you see it.

proteanplanet · 2023-07-12T14:37:05Z

I would like to suggest that rather than putting this into MPAS-Analysis, or perhaps in addition to putting this into MPAS-Analysis, we create an analysis member directly in MPAS-SeaIce that outputs the risk index.

xylar · 2023-07-15T15:58:32Z

@milenaveneziani, I'll be back from vacation on Tuesday and will take a look as soon as I can after that.

milenaveneziani · 2023-07-15T17:57:03Z

sounds good. Thanks @xylar.
There is one last thing to add to the title of the plot: a parameter value that comes from the config file, so that is super easy to do.

xylar

@torotorotaxi and @milenaveneziani, this looks promising!

I have a few suggested changes, most minor. The only more significant one would be to move the hard-coded arrays to a data file.

It would be nice if we could get rid of the rio>30 noisy results (darkest blue), but the things the we tried didn't work.

I think the problem is that things are either fractionally above or fractionally below 30.0, depending on machine roundoff. I would move the upper bound to either 30.01 or 29.99 so 30.0 is consistently on one side or the other.

xylar · 2023-07-18T16:00:49Z

mpas_analysis/sea_ice/climatology_map_risk_index_outcome.py

+# this are labels that should appear in the plot, to indicate the Polar Class of the vessel (included here for the moment)
+#       pic = ["PC1", "PC2", "PC3", "PC4", "PC5", "PC6", "PC7", "IA Super",\
+#              "IA", "IB", "IC", "Not Ice Strengthened"]
+
+        # reference floe thicknesses for calculation of Risk Index Values 
+        # (this values were agreed upon by Elizabeth Hunke, Andrew Roberts,
+        # and Gennaro D'Angelo based on literature and IMO description)
+        h_riv = np.array([0.5, 10, 15, 30, 50, 70, 100, 120, 170, 200, 250]) * 0.01 
+        # table of Risk Index Values (defined by IMO)
+        riv = np.array([[ 3, 3, 3, 3, 2, 2, 2, 2, 2, 2, 1, 1 ],\
+                        [ 3, 3, 3, 3, 2, 2, 2, 2, 2, 1, 1, 0 ],\
+                        [ 3, 3, 3, 3, 2, 2, 2, 2, 2, 1, 0,-1 ],\
+                        [ 3, 3, 3, 3, 2, 2, 2, 2, 1, 0,-1,-2 ],\
+                        [ 3, 3, 3, 3, 2, 2, 1, 1, 0,-1,-2,-2 ],\
+                        [ 3, 2, 2, 2, 2, 1, 1, 0,-1,-2,-3,-3 ],\
+                        [ 3, 2, 2, 2, 1, 1, 0,-1,-2,-3,-3,-3 ],\
+                        [ 3, 2, 2, 2, 2, 1, 0,-1,-2,-3,-4,-4 ],\
+                        [ 3, 2, 2, 2, 1, 0,-1,-2,-3,-4,-5,-5 ],\
+                        [ 3, 2, 2, 1, 0,-1,-2,-3,-4,-5,-6,-6 ],\
+                        [ 3, 2, 1, 0,-1,-2,-3,-4,-5,-6,-7,-8 ],\
+                        [ 3, 1, 0,-1,-2,-3,-4,-5,-6,-7,-8,-8 ]])
+


I'm pretty opposed to having arrays like this in the code, so it would be better to have all of this in a NetCDF or csv data file (with a date stamp). The file would also ideally be on the LCRC server and not part of the code. This would also allow us to update the data file with a new date stamp as needed.

xylar · 2023-07-18T17:39:34Z

mpas_analysis/sea_ice/climatology_map_risk_index_outcome.py

+        # sea-ice concentration conversion from range [0,1] to range [0,10]
+        scale_factor = 10
+        # polar class array index should be in the range [0,11], but not checked!
+        pc = self.config.getint(self.taskName, 'polarClass') - 1


It there a reason not to produce plots for all 12 classes? If yes, could there be a desire to produce plots for more than one class, in which case polarClass should become a list of polarClasses?

@xylar: should we put this file in a folder inside observations/SeaIce? If so, I need to figure out how to call it (will look at another task for reference) and remind myself how I update the documentation so that it's reflected in https://mpas-dev.github.io/MPAS-Analysis/develop/users_guide/observations.html#sea-ice-observations.

@milenaveneziani, please add information about your observations to the bottom of the sea-ice section of this XML file:
https://github.com/MPAS-Dev/MPAS-Analysis/blob/develop/mpas_analysis/obs/observational_datasets.xml#L1700
(before the icebergs section). You should be able to see what is expected from the other observations there.
Yes, please put your observations on LCRC under /lcrc/group/e3sm/public_html/diagnostics/observations/SeaIce in an appropriate subdirectory and, if it helps you with testing, you can make another copy under /lcrc/group/e3sm/diagnostics/observations/SeaIce. (Diagnostics get synced from public_html to that location as well as to other supported machines as part of E3SM-Unified testing and deployment.)

xylar · 2023-07-18T17:40:06Z

mpas_analysis/default.cfg

+# Polar Class of vessel according to IMO. Range is 1 to 12 (increments of 1).
+polarClass = 6


This should be documented a bit more to list the 12 polar classes, as is done in the code.

xylar · 2023-07-18T17:40:55Z

mpas_analysis/default.cfg

+# reference lat/lon for sea ice plots in the northern hemisphere
+minimumLatitude = -50
+referenceLongitude = 180


It seems like the comment need to be updated, since you have a minimum, rather than a reference, latitude.

Also, this is for the southern hemisphere.

xylar · 2023-07-18T17:41:36Z

mpas_analysis/default.cfg

+# Polar Class of vessel according to IMO. Range is 1 to 12 (increments of 1).
+polarClass = 6


This should be documented a bit more to list the 12 polar classes, as is done in the code.

xylar · 2023-07-18T17:41:53Z

mpas_analysis/default.cfg

+# reference lat/lon for sea ice plots in the northern hemisphere
+minimumLatitude = 50
+referenceLongitude = 0


It seems like the comment need to be updated, since you have a minimum, rather than a reference, latitude.

xylar · 2025-01-17T17:57:07Z

@milenaveneziani and @torotorotaxi, what would it take to get this finished and merge? it's been stale for about a year an a half so it's time to either push it through or let it go, I would say.

proteanplanet · 2025-01-17T19:04:15Z

My two cents worth - calculating the RIO as an MPAS-SeaIce analysis member has been on our list, and would provide a more accurate assessment of risk because it would access the sea ice thickness distribution in run time. I suggest dropping this PR and together pursuing the RIO output directly from MPAS-SeaIce instead.

milenaveneziani · 2025-01-17T20:50:56Z

@xylar, @proteanplanet: we actually need this for two projects with @stephenprice, so please do not close it. @torotorotaxi and I will discuss the best way forward for this PR shortly. Thanks.

stephenprice · 2025-01-30T22:38:31Z

@xylar -- Would it be possible for us to keep this PR open for a bit longer? We have an internal (LANL) project ongoing right now (winter-spring '25) that is working with/on this. Our assumption is that it will be in a position to merge following a bit more work during this time period. Thanks.

xylar · 2025-01-31T12:15:29Z

@stephenprice, I'm fine to leave this up. I was just checking a couple weeks back because it seemed abandoned.

milenaveneziani · 2025-02-04T16:36:44Z

I rebased and addressed a few minor comments about documentation in the config file. We are testing this now on chicoma

milenaveneziani · 2025-03-16T17:35:08Z

@xylar: for some reason Gennaro is unable to push to his own fork (permission issue..), which is weird, but in any case these are all his changes, which can be summarized as follows:

as per your request, the table with the RIV values has been moved to file. I decided to put it into the observations/SeaIce/RIO directory: let me know if you would prefer a different location.
the RIO calculation per se has been updated, adding an option to calculate RIO from the ice concentration categories. Plots have also been updated.

Things that for sure are left to do for this PR are: 1) add the RIV file to the public repo; 2) add something in the documentation that describes the RIV file.

milenaveneziani · 2025-03-16T17:37:33Z

One last thing that we need to decide with @stephenprice is whether we want to add something related to a specific route. I feel like there are two ways of doing it: 1) use the remapped RIO climatology that is calculated here and save the RIO-on-route information to file. This would require a simple addition to this PR; or 2) do it separately as a separate script using MPAS native fields (similarly to how I handle transect time series or vertical sections along a transect).
We can talk about this again after spring break.

milenaveneziani · 2025-03-16T17:40:41Z

Here are some sample plots. Using aggregated sea ice concentration:

milenaveneziani · 2025-03-16T17:43:01Z

and using the ice categories:

milenaveneziani · 2025-03-27T20:44:52Z

@xylar: an update on this. We have decided against adding anything else related to specific routes, so this PR should be almost ready to go (after resolving the conflicts).

The only things left are these:

add the RIV file to the public repo; 2) add something in the documentation that describes the RIV file.

Could you please remind me how to do the documentation update? (I will also take a look at the documentation for developers)

xylar · 2025-03-27T20:49:43Z

Here's the gist:
https://mpas-dev.github.io/MPAS-Analysis/develop/tutorials/dev_add_task.html#adding-documentation
Let me know if you have questions there.

xylar · 2025-03-27T20:52:06Z

For the RIV (which count as "obs", I think), I would appreciate an entry in https://github.com/MPAS-Dev/MPAS-Analysis/blob/develop/mpas_analysis/obs/observational_datasets.xml. Let me know if you don't think that's appropriate.

milenaveneziani · 2025-03-27T20:53:21Z

yes, that was my idea as well (to count them as obs)

milenaveneziani · 2025-03-31T22:46:42Z

@xylar: I followed the steps for adding documentation (plus I also added two list entries to docs/users_guide/analysis_tasks.rst, which was not mentioned in the developers documentation), but I don't think I can see the tasks listed in my test build of the docs. You'll probably know why right away..
I also added these files to the public repo:

pwd
/lcrc/group/e3sm/public_html/diagnostics/observations/SeaIce
chrlogin2:ac.milena/lcrc/.../SeaIce> ll RIV
total 3.5K
drwxrwsr-x 2 ac.milena E3SM 4.0K Mar 31 17:16 ./
drwxrwsr-x 9 ac.milena E3SM 4.0K Mar 31 15:28 ../
-rw-r-xr-- 1 ac.milena E3SM  633 Mar 31 17:13 obs.bib*
-rw-r-xr-- 1 ac.milena E3SM  874 Mar 31 17:16 README.md*
-rw-rw-r-- 1 ac.milena E3SM  681 Mar 31 15:30 riv_MSC.1_Circ.1519_6_June_2016.csv

xylar

@xylar: I followed the steps for adding documentation

Sounds good!

(plus I also added two list entries to docs/users_guide/analysis_tasks.rst, which was not mentioned in the developers documentation)

Sorry about that!

but I don't think I can see the tasks listed in my test build of the docs.

I build the docs locally and see the tasks fine. Maybe something went wrong in your build?

You'll probably know why right away..

I only see a few very minor issues, which I flagged below.

I also added these files to the public repo:

-rw-r-xr-- 1 ac.milena E3SM 633 Mar 31 17:13 obs.bib*
-rw-r-xr-- 1 ac.milena E3SM 874 Mar 31 17:16 README.md*

Could you chmod those two files to be group writable but not group executable? Not a big thing but presumably the intention.

xylar · 2025-04-01T02:31:28Z

docs/users_guide/tasks/climatologyMapRiskIndexOutcomeNH.rst

+climatologyMapRiskIndexOutcomeNH
+=======================================


Suggested change

climatologyMapRiskIndexOutcomeNH

=======================================

climatologyMapRiskIndexOutcomeNH

================================

xylar · 2025-04-01T02:31:44Z

docs/users_guide/tasks/climatologyMapRiskIndexOutcomeSH.rst

+climatologyMapRiskIndexOutcomeSH
+=======================================


Suggested change

climatologyMapRiskIndexOutcomeSH

=======================================

climatologyMapRiskIndexOutcomeSH

================================

xylar · 2025-04-01T02:32:26Z

docs/users_guide/tasks/climatologyMapRiskIndexOutcomeNH.rst

+Example Result
+--------------


Suggested change

Example Result

--------------

Observations

------------

:ref:`imo_riv`

Example Result

--------------

xylar · 2025-04-01T02:32:44Z

docs/users_guide/tasks/climatologyMapRiskIndexOutcomeSH.rst

+Example Result
+--------------


Suggested change

Example Result

--------------

Observations

------------

:ref:`imo_riv`

Example Result

--------------

milenaveneziani · 2025-04-01T03:46:29Z

I build the docs locally and see the tasks fine. Maybe something went wrong in your build?

it's possible. I'm glad you can see them just fine.

Could you chmod those two files to be group writable but not group executable?

yep, done.

I only see a few very minor issues, which I flagged below

I also just fixed these.

Not sure why there are conflicts in sea_ice/__init__.py: should I click on the web editor as suggested?

milenaveneziani · 2025-04-01T03:58:27Z

I think I have rebased successfully my local branch to origin/develop and solved the conflicts. I haven't pushed yet though, in case you prefer to handle things differently to preserve the commit history. Let me know what you think.

xylar · 2025-04-01T04:02:55Z

A rebase is my preferred way so go ahead.

milenaveneziani · 2025-04-01T13:22:12Z

This is ready to go on my side, @xylar. Let me know if there is anything else that should be done.

xylar · 2025-04-01T14:27:34Z

@milenaveneziani, I'm going to run my test suite and make sure I don't see any bugs.

One more question for you all. My reading is that this would be on by default and doesn't require any special analysis member or anything like that to be present in the simulation. I just want to make sure that we do want the new analysis on my default without users having to set any config options. (It wouldn't make it into the Unified release that might happen today or tomorrow, it's just too late for that, so it would start showing up in E3SM simulation results in 6 months or so.)

xylar · 2025-04-01T15:54:41Z

Here's what I'm seeing in the latest build of the docs. Looks good to me:
https://web.lcrc.anl.gov/public/e3sm/diagnostic_output/ac.xylar/analysis_testing/chrysalis/add_risk_index_outcome/docs/users_guide/tasks/climatologyMapRiskIndexOutcomeNH.html

milenaveneziani · 2025-04-01T16:02:43Z

My reading is that this would be on by default and doesn't require any special analysis member or anything like that to be present in the simulation

exactly. This only needs default sea ice output to run.

No problem about making it to e3sm-unified. We have our developer mpas conda builds that we can use for analysis.

And yes, the documentation also looks good.

xylar · 2025-04-01T16:59:31Z

I didn't think to sync the diagnostics fields before I ran my test suite, and my tests failed because the file is missing. It seems like others may run into a similar problem. Let me figure out why that's crashing MPAS-Analysis and not just causing the analysis task to fail.

xylar

@milenaveneziani and @torotorotaxi, I'm afraid I'm going to have to ask for some restructuring. It won't work to read the CSV file like you're doing in the constructor. The danger that this makes MPAS-Analysis as a whole unusable is just too high.

As an example, you can see how the Antarctic melt observations are read from a CSV file in run_task():
https://github.com/MPAS-Dev/MPAS-Analysis/blob/develop/mpas_analysis/ocean/climatology_map_antarctic_melt.py#L631

xylar · 2025-04-01T17:00:50Z

mpas_analysis/sea_ice/climatology_map_risk_index_outcome.py

+        polarclass = np.int_(polarclass) - 1
+
+        # read in table of Risk Index Values
+        #riv_csv = 'riv_MSC.1_Circ.1519_6_June_2016.csv'


Just a bit of suggested clean-up:

Suggested change

#riv_csv = 'riv_MSC.1_Circ.1519_6_June_2016.csv'

xylar · 2025-04-01T17:05:30Z

mpas_analysis/sea_ice/climatology_map_risk_index_outcome.py

+        IceClassLabels = np.genfromtxt(riv_csv, delimiter=',', skip_header=1,
+                                       dtype=str, usecols=(0,))


The problem here is that if riv_csv is missing, MPAS-Analysis will crash completely. That's not a good state to be in. This file needs to be read and used in a run_task() method, not in a constructor. As this is currently written, even if someone tries to disable this analysis, they are required to have this CSV file in order to run other analysis.

xylar · 2025-04-01T17:14:42Z

I can help with these changes if you need me to, but I'd be grateful if you could give it the first go.

milenaveneziani · 2025-04-01T17:45:57Z

hmm, not sure if I have the bandwidth to try this, definitely not this week.
Do you suppose a raise IOError could work in this case?

xylar · 2025-04-01T18:20:26Z

Do you suppose a raise IOError could work in this case?

That is what will happen if the file is missing anyway and what will crash all of MPAS-Analysis. The problem is that opening a diagnostics file that may or may not exist is simply not a safe thing to do in the constructor of an analysis task. You could have a try/except to catch this but then you would need to gracefully handle the case that the file doesn't exist. I think you would just ideally not be loading the file at all at construction time.

xylar · 2025-04-02T00:56:35Z

@milenaveneziani, I looked into it a bit more and there isn't a clean way to reorganize in the way I was hoping (so the file only gets read when the task gets run) but it's easy enough to just do nothing if the file is missing. I have some suggested clean up and a fix on my branch here:
https://github.com/xylar/MPAS-Analysis/tree/add_risk_index_outcome
Specifically, take a look at this addition:
xylar@e8355cd

Here's the clean-up I'd recommend:
xylar@947ba5d

If you're happy with both, please cherry-pick them (or hard-reset to my branch and force push).

xylar · 2025-04-02T02:24:17Z

We're getting there but the main vs. ctrl differences don't quite look right:

xylar · 2025-04-02T02:25:13Z

mpas_analysis/default.cfg

+# colormap for model/observations
+colormapNameResult = RdYlBu
+# whether the colormap is indexed or continuous
+colormapTypeResult = indexed
+# color indices into colormapName for filled contours
+colormapIndicesResult = [0, 56, 85, 170, 198, 227, 241, 248, 255, 255]
+# colormap levels/values for contour boundaries
+colorbarLevelsResult = numpy.linspace(-10., 30., 9)


You also need to define a colormap for the differences for main vs. control runs. (for NH and SH).

xylar · 2025-04-02T14:24:35Z

Here's an alternative approach I would recommend, in which the CSV file is stored in the python package instead of the diagnostics data. I realize I said above that the file would ideally be on the LCRC server. It seems like that approach isn't working out very well and that it needs to be in code to avoid issues when the data isn't found.

fcd818ada

xylar · 2025-04-03T14:19:07Z

I tested the alternative approach above. I had to make a small bug fix so I've updated the commit hash in the link. The results of my test suite with this approach look good:
https://web.lcrc.anl.gov/public/e3sm/diagnostic_output/ac.xasay-davis/analysis_testing/chrysalis/add_risk_index_outcome/

milenaveneziani · 2025-04-03T14:36:51Z

Sounds great!
So, I need to remove the files from the obs directory at LLRC and also the changes to the docs related to the obs part, right?

Will try to do this later this afternoon.
Thanks @xylar!

xylar · 2025-04-04T00:17:17Z

I think it would avoid confusing if we remove the LCRC RIV directory, yes. And, yes, please update the docs.

torotorotaxi changed the title ~~Add sea-ice risk index outcome~~ Add risk index outcome for navigation in ice-covered waters Jul 11, 2023

xylar reviewed Jul 18, 2023

View reviewed changes

milenaveneziani force-pushed the add_risk_index_outcome branch from c0c90e7 to 792484e Compare February 4, 2025 16:35

MPAS-Dev deleted a comment from torotorotaxi Feb 4, 2025

xylar reviewed Apr 1, 2025

View reviewed changes

“Gennaro” added 2 commits March 31, 2025 21:55

Add Risk Index Outcome analysis task

ceeebd5

Errors corrected and description added

a4d81ec

milenaveneziani added 3 commits March 31, 2025 21:56

Updated minor comments in config file

f312b83

Update of RIO computation. Moves RIV table to file.

cbfd36b

Adds RIO documentation and RIV entry in obs

5295bbb

milenaveneziani force-pushed the add_risk_index_outcome branch from 8b56fcb to 5295bbb Compare April 1, 2025 04:04

xylar requested changes Apr 1, 2025

View reviewed changes

xylar reviewed Apr 2, 2025

View reviewed changes

		# Polar Class of vessel according to IMO. Range is 1 to 12 (increments of 1).
		polarClass = 6

		climatologyMapRiskIndexOutcomeNH
		=======================================

		climatologyMapRiskIndexOutcomeSH
		=======================================

		IceClassLabels = np.genfromtxt(riv_csv, delimiter=',', skip_header=1,
		dtype=str, usecols=(0,))

Add risk index outcome for navigation in ice-covered waters #968

Are you sure you want to change the base?

Add risk index outcome for navigation in ice-covered waters #968

Uh oh!

Conversation

torotorotaxi commented Jul 11, 2023 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

torotorotaxi commented Jul 11, 2023

Uh oh!

torotorotaxi commented Jul 11, 2023

Uh oh!

torotorotaxi commented Jul 11, 2023

Uh oh!

milenaveneziani commented Jul 12, 2023

Uh oh!

proteanplanet commented Jul 12, 2023

Uh oh!

xylar commented Jul 15, 2023

Uh oh!

milenaveneziani commented Jul 15, 2023

Uh oh!

xylar left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

xylar commented Jan 17, 2025

Uh oh!

proteanplanet commented Jan 17, 2025

Uh oh!

milenaveneziani commented Jan 17, 2025

Uh oh!

stephenprice commented Jan 30, 2025

Uh oh!

xylar commented Jan 31, 2025

Uh oh!

milenaveneziani commented Feb 4, 2025

Uh oh!

milenaveneziani commented Mar 16, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

milenaveneziani commented Mar 16, 2025

Uh oh!

milenaveneziani commented Mar 16, 2025

Uh oh!

milenaveneziani commented Mar 16, 2025

Uh oh!

milenaveneziani commented Mar 27, 2025

Uh oh!

xylar commented Mar 27, 2025

Uh oh!

xylar commented Mar 27, 2025

Uh oh!

milenaveneziani commented Mar 27, 2025

Uh oh!

milenaveneziani commented Mar 31, 2025

Uh oh!

xylar left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

torotorotaxi commented Jul 11, 2023 •

edited

Loading

milenaveneziani commented Mar 16, 2025 •

edited

Loading

milenaveneziani commented Apr 1, 2025 •

edited

Loading

xylar Apr 1, 2025 •

edited

Loading

xylar commented Apr 2, 2025 •

edited

Loading