From bdda7ef55a3d63bc45330dc326d433eb9956a7e2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=E2=80=9CGennaro=E2=80=9D?= <“gennaro@lanl.gov”> Date: Fri, 23 Jun 2023 12:06:18 -0600 Subject: [PATCH 1/4] Add Risk Index Outcome analysis task --- mpas_analysis/__main__.py | 6 + mpas_analysis/default.cfg | 62 ++++ mpas_analysis/sea_ice/__init__.py | 2 + .../climatology_map_risk_index_outcome.py | 279 ++++++++++++++++++ 4 files changed, 349 insertions(+) create mode 100755 mpas_analysis/sea_ice/climatology_map_risk_index_outcome.py diff --git a/mpas_analysis/__main__.py b/mpas_analysis/__main__.py index d6bef5b5e..58ed3c219 100755 --- a/mpas_analysis/__main__.py +++ b/mpas_analysis/__main__.py @@ -251,6 +251,12 @@ def build_analysis_list(config, controlConfig): analyses.append(sea_ice.ClimatologyMapSeaIceProduction( config=config, mpas_climatology_task=seaIceClimatologyTask, hemisphere='SH', control_config=controlConfig)) + analyses.append(sea_ice.ClimatologyMapRiskIndexOutcoume( + config=config, mpas_climatology_task=seaIceClimatologyTask, + hemisphere='NH', control_config=controlConfig)) + analyses.append(sea_ice.ClimatologyMapRiskIndexOutcoume( + config=config, mpas_climatology_task=seaIceClimatologyTask, + hemisphere='SH', control_config=controlConfig)) analyses.append(sea_ice.ClimatologyMapSeaIceMelting( config=config, mpas_climatology_task=seaIceClimatologyTask, hemisphere='NH', control_config=controlConfig)) diff --git a/mpas_analysis/default.cfg b/mpas_analysis/default.cfg index 8173ccdd0..5f74931fa 100755 --- a/mpas_analysis/default.cfg +++ b/mpas_analysis/default.cfg @@ -3964,6 +3964,68 @@ volNH = PIOMAS/PIOMASvolume_monthly_climo_20180710.nc volSH = none +[f'climatologyMapRiskIndexOutcomeNH] +## options related to plotting maps of risk index outcomes for navigability +## in sea-ice covered water from climatologies + +# colormap for model/observations +colormapNameResult = RdYlBu +# whether the colormap is indexed or continuous +colormapTypeResult = indexed +# color indices into colormapName for filled contours +colormapIndicesResult = numpy.array(numpy.linspace(0, 255, 9), int) +# colormap levels/values for contour boundaries +colorbarLevelsResult = numpy.linspace(-10., 30., 8) + +# Months or seasons to plot (Jan, Feb, Mar, Apr, May, Jun, Jul, Aug, Sep, Oct, +# Nov, Dec, JFM, AMJ, JAS, OND, ANN) +seasons = ['MAM','SON'] + +# comparison grid(s) on which to plot analysis +comparisonGrids = ['arctic'] + +# Polar Class of vessel according to IMO. Range is 1 to 12 (increments of 1). +polarClass = 6 + +# reference lat/lon for sea ice plots in the northern hemisphere +minimumLatitude = 50 +referenceLongitude = 0 + +# arrange subplots vertically? +#vertical = False + + +[f'climatologyMapRiskIndexOutcomeSH] +## options related to plotting maps of risk index outcomes for navigability +## in sea-ice covered water from climatologies + +# colormap for model/observations +colormapNameResult = RdYlBu +# whether the colormap is indexed or continuous +colormapTypeResult = indexed +# color indices into colormapName for filled contours +colormapIndicesResult = numpy.array(numpy.linspace(0, 255, 9), int) +# colormap levels/values for contour boundaries +colorbarLevelsResult = numpy.linspace(-10., 30., 8) + +# Months or seasons to plot (Jan, Feb, Mar, Apr, May, Jun, Jul, Aug, Sep, Oct, +# Nov, Dec, JFM, AMJ, JAS, OND, ANN) +seasons = ['MAM','SON'] + +# comparison grid(s) on which to plot analysis +comparisonGrids = ['antarctic'] + +# Polar Class of vessel according to IMO. Range is 1 to 12 (increments of 1). +polarClass = 6 + +# reference lat/lon for sea ice plots in the northern hemisphere +minimumLatitude = -50 +referenceLongitude = 180 + +# arrange subplots vertically? +#vertical = False + + [climatologyMapSeaIceProductionNH] # options related to plotting horizontally remapped climatologies of # sea ice production against control model results and observations diff --git a/mpas_analysis/sea_ice/__init__.py b/mpas_analysis/sea_ice/__init__.py index 300523152..430865414 100644 --- a/mpas_analysis/sea_ice/__init__.py +++ b/mpas_analysis/sea_ice/__init__.py @@ -9,3 +9,5 @@ ClimatologyMapSeaIceProduction from mpas_analysis.sea_ice.climatology_map_melting import \ ClimatologyMapSeaIceMelting +from mpas_analysis.sea_ice.climatology_map_risk_index_outcome import \ + ClimatologyMapRiskIndexOutcoume diff --git a/mpas_analysis/sea_ice/climatology_map_risk_index_outcome.py b/mpas_analysis/sea_ice/climatology_map_risk_index_outcome.py new file mode 100755 index 000000000..16c0699cc --- /dev/null +++ b/mpas_analysis/sea_ice/climatology_map_risk_index_outcome.py @@ -0,0 +1,279 @@ +# Copyright (c) 2017, Los Alamos National Security, LLC (LANS) +# and the University Corporation for Atmospheric Research (UCAR). +# +# Unless noted otherwise source code is licensed under the BSD license. +# Additional copyright and license information can be found in the LICENSE file +# distributed with this code, or at http://mpas-dev.github.com/license.html +# + +import xarray as xr +from pyremap import LatLon2DGridDescriptor + +from mpas_analysis.shared import AnalysisTask + +from mpas_analysis.shared.climatology import RemapMpasClimatologySubtask, \ + RemapObservedClimatologySubtask + +from mpas_analysis.shared.plot import PlotClimatologyMapSubtask + +from mpas_analysis.shared.io.utility import build_obs_path + + +class ClimatologyMapRiskIndexOutcoume(AnalysisTask): + """ + An analysis task for evaluating the Risk Index Outcome + for navigation in sea-ice covered water. + (https://www.imorules.com/GUID-2C1D86CB-5D58-490F-B4D4-46C057E1D102.html) + """ + # Authors + # ------- + # Gennaro D'Angelo, Milena Veneziani + + def __init__(self, config, mpas_climatology_task, hemisphere, + control_config=None): + """ + Construct the analysis task. + + Parameters + ---------- + config : mpas_tools.config.MpasConfigParser + Configuration options + + mpas_climatology_task : mpas_analysis.shared.climatology.MpasClimatologyTask + The task that produced the climatology to be remapped and plotted + + hemisphere : {'NH', 'SH'} + The hemisphere to plot + + control_config : mpas_tools.config.MpasConfigParser, optional + Configuration options for a control run (if any) + """ + # Authors + # ------- + # Gennaro D'Angelo, Milena Veneziani + + task_name = f'climatologyMapRiskIndexOutcome{hemisphere}' + + field_name = 'riskIndexOutcome' + + tags = ['climatology', 'horizontalMap', field_name] + if hemisphere == 'NH': + tags = tags + ['arctic'] + else: + tags = tags + ['antarctic'] + + # call the constructor from the base class (AnalysisTask) + super().__init__(config=config, taskName=task_name, + componentName='seaIce', tags=tags) + + self.mpas_climatology_task = mpas_climatology_task + + section_name = self.taskName + + if hemisphere == 'NH': + hemisphere_long= 'Northern' + else: + hemisphere_long= 'Southern' + + # read in what seasons we want to plot + seasons = config.getexpression(section_name, 'seasons') + + if len(seasons) == 0: + raise ValueError(f'config section {section_name} does not contain ' + f'valid list of seasons') + + comparison_grid_names = config.getexpression(section_name, + 'comparisonGrids') + + if len(comparison_grid_names) == 0: + raise ValueError(f'config section {section_name} does not contain ' + f'valid list of comparison grids') + + variable_list = ['timeMonthly_avg_iceAreaCell', + 'timeMonthly_avg_iceVolumeCell'] + + remap_climatology_subtask = RemapMpasRiskIndexOutcomeClimatology( + mpas_climatology_task=mpas_climatology_task, + parent_task=self, + climatology_name=f'{field_name}{hemisphere}', + variable_list=variable_list, + comparison_grid_names=comparison_grid_names, + seasons=seasons) + + self.add_subtask(remap_climatology_subtask) + + for season in seasons: + for comparison_grid_name in comparison_grid_names: + + if control_config is None: + remap_observations_subtask = None + gallery_name = None + ref_title_label = None + ref_field_name = None + diff_title_label = 'Model - Observations' + + else: + control_run_name = control_config.get('runs', 'mainRunName') + gallery_name = None + ref_title_label = f'Control: {control_run_name}' + field_name = field_name + diff_title_label = 'Main - Control' + + image_caption = f'{season} Climatology Map of ' \ + f'{hemisphere_long}-Hemisphere Risk Index Outcome' + gallery_group = f'{hemisphere_long}-Hemisphere Risk Index Outcome' + # make a new subtask for this season and comparison grid + subtask = PlotClimatologyMapSubtask( + parentTask=self, season=season, + comparisonGridName=comparison_grid_name, + remapMpasClimatologySubtask=remap_climatology_subtask, + remapObsClimatologySubtask=None, + controlConfig=control_config) + + subtask.set_plot_info( + outFileLabel=f'risk_index_outcome{hemisphere}', + fieldNameInTitle='Risk Index Outcome', + mpasFieldName=field_name, + refFieldName=field_name, + refTitleLabel=ref_title_label, + diffTitleLabel=diff_title_label, + unitsLabel=r'', + imageCaption=image_caption, + galleryGroup=gallery_group, + groupSubtitle=None, + groupLink=f'{hemisphere.lower()}_rio', + galleryName=gallery_name, + extend=None) + + self.add_subtask(subtask) + + +class RemapMpasRiskIndexOutcomeClimatology(RemapMpasClimatologySubtask): + """ + A subtask for computing climatologies of risk index outcome from + climatologies of sea-ice concentration and thickness. + """ + def __init__(self, mpas_climatology_task, parent_task, climatology_name, + variable_list, seasons, comparison_grid_names): + + """ + Construct the analysis task and adds it as a subtask of the + ``parent_task``. + Parameters + ---------- + mpas_climatology_task : mpas_analysis.shared.climatology.MpasClimatologyTask + The task that produced the climatology to be remapped + parent_task : mpas_analysis.shared.AnalysisTask + The parent task, used to get the ``taskName``, ``config`` and + ``componentName`` + climatology_name : str + A name that describes the climatology (e.g. a short version of + the important field(s) in the climatology) used to name the + subdirectories for each stage of the climatology + variable_list : list of str + A list of variable names in ``timeSeriesStatsMonthly`` to be + included in the climatologies + seasons : list of str, optional + A list of seasons (keys in ``shared.constants.monthDictionary``) + to be computed or ['none'] (not ``None``) if only monthly + climatologies are needed. + comparison_grid_names : list of {'latlon', 'antarctic'} + The name(s) of the comparison grid to use for remapping. + """ + + subtask_name = f'remapMpasClimatology_RiskIndexOutcome' + # call the constructor from the base class + # (RemapMpasClimatologySubtask) + super().__init__( + mpas_climatology_task, parent_task, climatology_name, + variable_list, seasons, comparison_grid_names) + + self.mpas_climatology_task = mpas_climatology_task + self.variable_list = variable_list + + def setup_and_check(self): + """ + Perform steps to set up the analysis and check for errors in the setup. + """ + + # first, call setup_and_check from the base class + # (RemapMpasClimatologySubtask), which will set up remappers and add + # variables to mpas_climatology_task + super().setup_and_check() + + # don't add the variables and seasons to mpas_climatology_task until + # we're sure this subtask is supposed to run + self.mpas_climatology_task.add_variables(self.variable_list, + self.seasons) + + def customize_masked_climatology(self, climatology, season): + """ + Compute the risk index outcome from sea-ice concentration and thickness. + + Parameters + ---------- + climatology : xarray.Dataset + the climatology data set + season : str + The name of the season to be masked + Returns + ------- + climatology : xarray.Dataset + the modified climatology data set + """ + + rio = self._compute_risk_index_outcome(climatology) + + climatology['riskIndexOutcome'] = rio + climatology.riskIndexOutcome.attrs['units'] = '' + climatology = climatology.drop_vars(self.variable_list) + + return climatology + + def _compute_risk_index_outcome(self, climatology): + """ + Compute the risk index outcome from sea-ice concentration and thickness. + (https://www.imorules.com/GUID-2C1D86CB-5D58-490F-B4D4-46C057E1D102.html) + """ + ds_restart = xr.open_dataset(self.restartFileName) + ds_restart = ds_restart.isel(Time=0) + + scale_factor = 10 + pc = self.config.get(self.taskName, 'polarClass') - 1 + + concentration = climatology['timeMonthly_avg_iceAreaCell'] + thickness = climatology['timeMonthly_avg_iceVolumeCell'] + +# pic = ["PC1", "PC2", "PC3", "PC4", "PC5", "PC6", "PC7", "IA Super",\ +# "IA", "IB", "IC", "Not Ice Strengthened"] + + h_riv = np.array([0.5, 10, 15, 30, 50, 70, 100, 120, 170, 200, 250]) * 0.01 + 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 ]]) + + riv_iceCell = np.nan*np.ones(np.shape(thickness)) + + riv_mask = np.where(thickness < h_riv[0]) + riv_iceCell[riv_mask] = riv[pc, 0] + + for ind in range(len(h_riv)-1): +# riv_mask = np.logical_and(thickness >= h_riv[ind], thickness < h_riv[ind+1]) + riv_mask = np.where(thickness >= h_riv[ind]) + riv_iceCell[riv_mask] = riv[pc, ind+1] + + riv_mask = np.where(thickness >= h_riv[-1]) + riv_iceCell[riv_mask] = riv[pc, -1] + + rio = ((1.0 - concentration) * riv[pc, 0] + concentration * riv_iceCell) * units_scale_factor + + return rio From c288459716704d92ec88d440846d22f9396bde9c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=E2=80=9CGennaro=E2=80=9D?= <“gennaro@lanl.gov”> Date: Fri, 23 Jun 2023 16:54:50 -0600 Subject: [PATCH 2/4] Typos Corrected --- mpas_analysis/__main__.py | 4 ++-- mpas_analysis/default.cfg | 4 ++-- mpas_analysis/sea_ice/__init__.py | 2 +- mpas_analysis/sea_ice/climatology_map_risk_index_outcome.py | 2 +- 4 files changed, 6 insertions(+), 6 deletions(-) diff --git a/mpas_analysis/__main__.py b/mpas_analysis/__main__.py index 58ed3c219..9a92da749 100755 --- a/mpas_analysis/__main__.py +++ b/mpas_analysis/__main__.py @@ -251,10 +251,10 @@ def build_analysis_list(config, controlConfig): analyses.append(sea_ice.ClimatologyMapSeaIceProduction( config=config, mpas_climatology_task=seaIceClimatologyTask, hemisphere='SH', control_config=controlConfig)) - analyses.append(sea_ice.ClimatologyMapRiskIndexOutcoume( + analyses.append(sea_ice.ClimatologyMapRiskIndexOutcome( config=config, mpas_climatology_task=seaIceClimatologyTask, hemisphere='NH', control_config=controlConfig)) - analyses.append(sea_ice.ClimatologyMapRiskIndexOutcoume( + analyses.append(sea_ice.ClimatologyMapRiskIndexOutcome( config=config, mpas_climatology_task=seaIceClimatologyTask, hemisphere='SH', control_config=controlConfig)) analyses.append(sea_ice.ClimatologyMapSeaIceMelting( diff --git a/mpas_analysis/default.cfg b/mpas_analysis/default.cfg index 5f74931fa..bcd050931 100755 --- a/mpas_analysis/default.cfg +++ b/mpas_analysis/default.cfg @@ -3964,7 +3964,7 @@ volNH = PIOMAS/PIOMASvolume_monthly_climo_20180710.nc volSH = none -[f'climatologyMapRiskIndexOutcomeNH] +[climatologyMapRiskIndexOutcomeNH] ## options related to plotting maps of risk index outcomes for navigability ## in sea-ice covered water from climatologies @@ -3995,7 +3995,7 @@ referenceLongitude = 0 #vertical = False -[f'climatologyMapRiskIndexOutcomeSH] +[climatologyMapRiskIndexOutcomeSH] ## options related to plotting maps of risk index outcomes for navigability ## in sea-ice covered water from climatologies diff --git a/mpas_analysis/sea_ice/__init__.py b/mpas_analysis/sea_ice/__init__.py index 430865414..cbe635469 100644 --- a/mpas_analysis/sea_ice/__init__.py +++ b/mpas_analysis/sea_ice/__init__.py @@ -10,4 +10,4 @@ from mpas_analysis.sea_ice.climatology_map_melting import \ ClimatologyMapSeaIceMelting from mpas_analysis.sea_ice.climatology_map_risk_index_outcome import \ - ClimatologyMapRiskIndexOutcoume + ClimatologyMapRiskIndexOutcome diff --git a/mpas_analysis/sea_ice/climatology_map_risk_index_outcome.py b/mpas_analysis/sea_ice/climatology_map_risk_index_outcome.py index 16c0699cc..0aaba2cf1 100755 --- a/mpas_analysis/sea_ice/climatology_map_risk_index_outcome.py +++ b/mpas_analysis/sea_ice/climatology_map_risk_index_outcome.py @@ -19,7 +19,7 @@ from mpas_analysis.shared.io.utility import build_obs_path -class ClimatologyMapRiskIndexOutcoume(AnalysisTask): +class ClimatologyMapRiskIndexOutcome(AnalysisTask): """ An analysis task for evaluating the Risk Index Outcome for navigation in sea-ice covered water. From 502fda8c5f063d540081fd9f042f37a179a6a6ad Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=E2=80=9CGennaro=E2=80=9D?= <“gennaro@lanl.gov”> Date: Mon, 26 Jun 2023 11:57:00 -0600 Subject: [PATCH 3/4] Errors corrected and description added --- mpas_analysis/default.cfg | 12 +++--- .../climatology_map_risk_index_outcome.py | 37 +++++++++++++++---- 2 files changed, 35 insertions(+), 14 deletions(-) diff --git a/mpas_analysis/default.cfg b/mpas_analysis/default.cfg index bcd050931..2a35cd571 100755 --- a/mpas_analysis/default.cfg +++ b/mpas_analysis/default.cfg @@ -3973,13 +3973,13 @@ colormapNameResult = RdYlBu # whether the colormap is indexed or continuous colormapTypeResult = indexed # color indices into colormapName for filled contours -colormapIndicesResult = numpy.array(numpy.linspace(0, 255, 9), int) +colormapIndicesResult = numpy.array(numpy.linspace(0, 255, 10), int) # colormap levels/values for contour boundaries -colorbarLevelsResult = numpy.linspace(-10., 30., 8) +colorbarLevelsResult = numpy.linspace(-10., 30., 9) # Months or seasons to plot (Jan, Feb, Mar, Apr, May, Jun, Jul, Aug, Sep, Oct, # Nov, Dec, JFM, AMJ, JAS, OND, ANN) -seasons = ['MAM','SON'] +seasons = ['AMJ','OND'] # comparison grid(s) on which to plot analysis comparisonGrids = ['arctic'] @@ -4004,13 +4004,13 @@ colormapNameResult = RdYlBu # whether the colormap is indexed or continuous colormapTypeResult = indexed # color indices into colormapName for filled contours -colormapIndicesResult = numpy.array(numpy.linspace(0, 255, 9), int) +colormapIndicesResult = numpy.array(numpy.linspace(0, 255, 10), int) # colormap levels/values for contour boundaries -colorbarLevelsResult = numpy.linspace(-10., 30., 8) +colorbarLevelsResult = numpy.linspace(-10., 30., 9) # Months or seasons to plot (Jan, Feb, Mar, Apr, May, Jun, Jul, Aug, Sep, Oct, # Nov, Dec, JFM, AMJ, JAS, OND, ANN) -seasons = ['MAM','SON'] +seasons = ['AMJ','OND'] # comparison grid(s) on which to plot analysis comparisonGrids = ['antarctic'] diff --git a/mpas_analysis/sea_ice/climatology_map_risk_index_outcome.py b/mpas_analysis/sea_ice/climatology_map_risk_index_outcome.py index 0aaba2cf1..cc656fd60 100755 --- a/mpas_analysis/sea_ice/climatology_map_risk_index_outcome.py +++ b/mpas_analysis/sea_ice/climatology_map_risk_index_outcome.py @@ -7,6 +7,7 @@ # import xarray as xr +import numpy as np from pyremap import LatLon2DGridDescriptor from mpas_analysis.shared import AnalysisTask @@ -208,7 +209,7 @@ def setup_and_check(self): def customize_masked_climatology(self, climatology, season): """ - Compute the risk index outcome from sea-ice concentration and thickness. + Compute the risk index outcome from sea-ice concentration and (floe) thickness. Parameters ---------- @@ -232,22 +233,27 @@ def customize_masked_climatology(self, climatology, season): def _compute_risk_index_outcome(self, climatology): """ - Compute the risk index outcome from sea-ice concentration and thickness. + Compute the risk index outcome from sea-ice concentration and (floe) thickness, + as outlined in by International Maritime Organization (IMO) document. (https://www.imorules.com/GUID-2C1D86CB-5D58-490F-B4D4-46C057E1D102.html) """ ds_restart = xr.open_dataset(self.restartFileName) ds_restart = ds_restart.isel(Time=0) + # sea-ice concentration conversion from range [0,1] to range [0,10] scale_factor = 10 - pc = self.config.get(self.taskName, 'polarClass') - 1 - - concentration = climatology['timeMonthly_avg_iceAreaCell'] - thickness = climatology['timeMonthly_avg_iceVolumeCell'] + # polar class array index should be in the range [0,11], but not checked! + pc = self.config.getint(self.taskName, 'polarClass') - 1 +# 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 ],\ @@ -261,19 +267,34 @@ def _compute_risk_index_outcome(self, climatology): [ 3, 2, 1, 0,-1,-2,-3,-4,-5,-6,-7,-8 ],\ [ 3, 1, 0,-1,-2,-3,-4,-5,-6,-7,-8,-8 ]]) - riv_iceCell = np.nan*np.ones(np.shape(thickness)) + concentration = climatology['timeMonthly_avg_iceAreaCell'] + thickness = climatology['timeMonthly_avg_iceVolumeCell'] + # these lines may not be required, but rio should not be used + # to check concentration and thickness + concentration[np.where(concentration<0.0)] = 0.0 + concentration[np.where(concentration>1.0)] = 1.0 + thickness[np.where(concentration==0.0)] = 0.0 + + # introduce riv array and set values in open water + riv_iceCell = np.nan*np.ones(np.shape(thickness)) riv_mask = np.where(thickness < h_riv[0]) riv_iceCell[riv_mask] = riv[pc, 0] + # set riv values for chosen Polar Class of vessel for ind in range(len(h_riv)-1): # riv_mask = np.logical_and(thickness >= h_riv[ind], thickness < h_riv[ind+1]) riv_mask = np.where(thickness >= h_riv[ind]) riv_iceCell[riv_mask] = riv[pc, ind+1] + # set riv for highest floe thickness riv_mask = np.where(thickness >= h_riv[-1]) riv_iceCell[riv_mask] = riv[pc, -1] - rio = ((1.0 - concentration) * riv[pc, 0] + concentration * riv_iceCell) * units_scale_factor + # Risk Index Outcome for single-category ice. There are only + # two terms per cell: one coming from the fraction of the cell + # occupied by open water and one coming from the fraction occupied + # by sea ice (rio <= 30 by IMO definition) + rio = scale_factor * ((1.0 - concentration) * riv[pc, 0] + concentration * riv_iceCell) return rio From c0c90e7ba81eb15ebb2f0c3fe8c8a055aa3c8250 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=E2=80=9CGennaro=E2=80=9D?= <“gennaro@lanl.gov”> Date: Mon, 26 Jun 2023 12:07:23 -0600 Subject: [PATCH 4/4] Some typos corrected --- mpas_analysis/sea_ice/climatology_map_risk_index_outcome.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/mpas_analysis/sea_ice/climatology_map_risk_index_outcome.py b/mpas_analysis/sea_ice/climatology_map_risk_index_outcome.py index cc656fd60..e739f9756 100755 --- a/mpas_analysis/sea_ice/climatology_map_risk_index_outcome.py +++ b/mpas_analysis/sea_ice/climatology_map_risk_index_outcome.py @@ -234,7 +234,7 @@ def customize_masked_climatology(self, climatology, season): def _compute_risk_index_outcome(self, climatology): """ Compute the risk index outcome from sea-ice concentration and (floe) thickness, - as outlined in by International Maritime Organization (IMO) document. + as outlined in the International Maritime Organization (IMO) document. (https://www.imorules.com/GUID-2C1D86CB-5D58-490F-B4D4-46C057E1D102.html) """ ds_restart = xr.open_dataset(self.restartFileName) @@ -251,7 +251,7 @@ def _compute_risk_index_outcome(self, climatology): # 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) + # 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 ],\ @@ -293,7 +293,7 @@ def _compute_risk_index_outcome(self, climatology): # Risk Index Outcome for single-category ice. There are only # two terms per cell: one coming from the fraction of the cell - # occupied by open water and one coming from the fraction occupied + # covered by open water and one coming from the fraction covered # by sea ice (rio <= 30 by IMO definition) rio = scale_factor * ((1.0 - concentration) * riv[pc, 0] + concentration * riv_iceCell)