-
Notifications
You must be signed in to change notification settings - Fork 52
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
Add risk index outcome for navigation in ice-covered waters #968
base: develop
Are you sure you want to change the base?
Changes from all commits
File filter
Filter by extension
Conversations
Jump to
Diff view
Diff view
There are no files selected for viewing
Original file line number | Diff line number | Diff line change |
---|---|---|
|
@@ -3964,6 +3964,68 @@ volNH = PIOMAS/PIOMASvolume_monthly_climo_20180710.nc | |
volSH = none | ||
|
||
|
||
[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, 10), int) | ||
# colormap levels/values for contour boundaries | ||
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 = ['AMJ','OND'] | ||
|
||
# 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 | ||
Comment on lines
+3990
to
+3992
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. It seems like the comment need to be updated, since you have a minimum, rather than a reference, latitude. |
||
|
||
# arrange subplots vertically? | ||
#vertical = False | ||
|
||
|
||
[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, 10), int) | ||
# colormap levels/values for contour boundaries | ||
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 = ['AMJ','OND'] | ||
|
||
# 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 | ||
Comment on lines
+4018
to
+4019
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. This should be documented a bit more to list the 12 polar classes, as is done in the code. |
||
|
||
# reference lat/lon for sea ice plots in the northern hemisphere | ||
minimumLatitude = -50 | ||
referenceLongitude = 180 | ||
Comment on lines
+4021
to
+4023
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. It seems like the comment need to be updated, since you have a minimum, rather than a reference, latitude. There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Also, this is for the southern hemisphere. |
||
|
||
# arrange subplots vertically? | ||
#vertical = False | ||
|
||
|
||
[climatologyMapSeaIceProductionNH] | ||
# options related to plotting horizontally remapped climatologies of | ||
# sea ice production against control model results and observations | ||
|
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,300 @@ | ||
# 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 | ||
import numpy as np | ||
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 ClimatologyMapRiskIndexOutcome(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 (floe) 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 (floe) thickness, | ||
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) | ||
ds_restart = ds_restart.isel(Time=0) | ||
|
||
# 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 | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. 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 |
||
|
||
# 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 ]]) | ||
|
||
Comment on lines
+248
to
+269
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. 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. |
||
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] | ||
|
||
# Risk Index Outcome for single-category ice. There are only | ||
# two terms per cell: one coming from the fraction of the cell | ||
# 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) | ||
|
||
return rio |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This should be documented a bit more to list the 12 polar classes, as is done in the code.