Update to compass 1.4.0-alpha.4, numpy >=2.0, mpas_tools 0.34.1, matplotlib >=3.9.0

compass/landice/tests/calving_dt_convergence/dt_convergence_test/__init__.py

@@ -1,10 +1,11 @@
-from compass.testcase import TestCase
-from compass.landice.tests.calving_dt_convergence.run_model import RunModel
+import matplotlib.cm
+import matplotlib.pyplot as plt
+import netCDF4
 # from compass.validate import compare_variables  # not currently used
 import numpy
-import netCDF4
-import matplotlib.pyplot as plt
-import matplotlib.cm
+
+from compass.landice.tests.calving_dt_convergence.run_model import RunModel
+from compass.testcase import TestCase
 
 
 class DtConvergenceTest(TestCase):
@@ -28,7 +29,7 @@ def __init__(self, test_group, mesh, calving, velo):
         test_group : compass.landice.tests.calving_dt_convergence.CalvingDtConvergence
             The test group that this test case belongs to
 
-        """
+        """  # noqa: E501
         self.name = f'calving_dt_convergence_test_{mesh}_{calving}_{velo}'
         subdir = f'{mesh}.{calving}.{velo}'
         super().__init__(test_group=test_group, name=self.name, subdir=subdir)
@@ -77,7 +78,8 @@ def validate(self):
         ax[3].set_ylabel('# warnings', color='c')
         ax2 = ax[3].twinx()
         ax2.set_ylabel('fraction with warnings', color='g')
-        colors = matplotlib.cm.jet(numpy.linspace(0, 1, len(self.fractions)))
+        jet = matplotlib.colormaps['jet']
+        colors = jet(numpy.linspace(0, 1, len(self.fractions)))
         nWarn = numpy.zeros([len(self.fractions)])
         nTimesteps = numpy.zeros([len(self.fractions)])
 
@@ -91,7 +93,7 @@ def validate(self):
             ax[0].plot(yr[1:], calv[1:], '-', label=f'{frac:.2f}',
                        color=colors[i])
 
-            ax[1].plot(yr[1:], (calv[1:]*deltat[1:]).cumsum(), '-',
+            ax[1].plot(yr[1:], (calv[1:] * deltat[1:]).cumsum(), '-',
                        color=colors[i])
 
             ratio = f.variables['dtCalvingCFLratio'][:]
@@ -107,7 +109,7 @@ def validate(self):
             nWarn[i] = logcontents.count("WARNING: Failed to ablate")
             nTimesteps[i] = logcontents.count("Starting timestep number")
             ax[3].plot(frac, nWarn[i], 'co')
-            ax2.plot(frac, nWarn[i]/nTimesteps[i], 'gx')
+            ax2.plot(frac, nWarn[i] / nTimesteps[i], 'gx')
 
             f.close()
             i += 1

compass/landice/tests/eismint2/standard_experiments/visualize.py

@@ -1,6 +1,7 @@
 import datetime
-import netCDF4
+
 import matplotlib.pyplot as plt
+import netCDF4
 import numpy as np
 from scipy.interpolate import griddata
 
@@ -19,7 +20,7 @@ def __init__(self, test_case):
         ----------
         test_case : compass.landice.tests.eismint2.standard_experiments.StandardExperiments
             The test case this step belongs to
-        """
+        """  # noqa: E501
         super().__init__(test_case=test_case, name='visualize')
 
         # depending on settings, this may produce no outputs, so we won't add
@@ -70,8 +71,8 @@ def visualize_eismint2(config, logger, experiment):
 
     # open supplied MPAS output file and get variables needed
     filein = netCDF4.Dataset(filename, 'r')
-    xCell = filein.variables['xCell'][:]/1000.0
-    yCell = filein.variables['yCell'][:]/1000.0
+    xCell = filein.variables['xCell'][:] / 1000.0
+    yCell = filein.variables['yCell'][:] / 1000.0
     xtime = filein.variables['xtime'][:]
     nCells = len(filein.dimensions['nCells'])
     nVertLevels = len(filein.dimensions['nVertLevels'])
@@ -103,16 +104,17 @@ def visualize_eismint2(config, logger, experiment):
 
     # make an educated guess about how big the markers should be.
     if nCells**0.5 < 100.0:
-        markersize = max(int(round(3600.0/(nCells**0.5))), 1)
+        markersize = max(int(round(3600.0 / (nCells**0.5))), 1)
         # use hexes if the points are big enough, otherwise just dots
         markershape = 'h'
     else:
-        markersize = max(int(round(1800.0/(nCells**0.5))), 1)
+        markersize = max(int(round(1800.0 / (nCells**0.5))), 1)
         markershape = '.'
     logger.info('Using a markersize of {}'.format(markersize))
 
     fig = plt.figure(1, facecolor='w')
-    fig.suptitle('Payne et al. Fig. 1, 3, 6, 9, or 11', fontsize=10, fontweight='bold')
+    fig.suptitle('Payne et al. Fig. 1, 3, 6, 9, or 11', fontsize=10,
+                 fontweight='bold')
 
     iceIndices = np.where(thickness[timelev, :] > 10.0)[0]
     plt.scatter(xCell[iceIndices], yCell[iceIndices], markersize,
@@ -138,7 +140,8 @@ def visualize_eismint2(config, logger, experiment):
         plt.savefig('EISMINT2-{}-basaltemp.png'.format(experiment), dpi=150)
 
     # ================
-    # STEADY STATE MAPS -  panels b and c are switched and with incorrect units in the paper
+    # STEADY STATE MAPS -  panels b and c are switched and with incorrect
+    #                      units in the paper
     # ================
     fig = plt.figure(2, facecolor='w', figsize=(12, 6), dpi=72)
     fig.suptitle('Payne et al. Fig. 2 or 4', fontsize=10, fontweight='bold')
@@ -152,7 +155,7 @@ def visualize_eismint2(config, logger, experiment):
                 edgecolors='none')
 
     # add contours of ice thickness over the top
-    contour_intervals = np.linspace(0.0, 5000.0,  int(5000.0/250.0)+1)
+    contour_intervals = np.linspace(0.0, 5000.0, int(5000.0 / 250.0) + 1)
     _contour_mpas(thickness[timelev, :], nCells, xCell, yCell,
                   contour_levs=contour_intervals)
 
@@ -167,17 +170,17 @@ def visualize_eismint2(config, logger, experiment):
 
     flux = np.zeros((nCells,))
     for k in range(nVertLevels):
-        speedLevel = (uReconstructX[timelev, :, k:k+2].mean(axis=1)**2 +
-                      uReconstructY[timelev, :, k:k+2].mean(axis=1)**2)**0.5
+        speedLevel = (uReconstructX[timelev, :, k:k + 2].mean(axis=1)**2 +
+                      uReconstructY[timelev, :, k:k + 2].mean(axis=1)**2)**0.5
         flux += speedLevel * thickness[timelev, :] * layerThicknessFractions[k]
 
     plt.scatter(xCell[iceIndices], yCell[iceIndices], markersize,
                 c=np.array([[0.8, 0.8, 0.8], ]), marker=markershape,
                 edgecolors='none')
 
     # add contours over the top
-    contour_intervals = np.linspace(0.0, 20.0,  11)
-    _contour_mpas(flux * 3600.0*24.0*365.0 / 10000.0, nCells, xCell, yCell,
+    contour_intervals = np.linspace(0.0, 20.0, 11)
+    _contour_mpas(flux * 3600.0 * 24.0 * 365.0 / 10000.0, nCells, xCell, yCell,
                   contour_levs=contour_intervals)
     ax.set_aspect('equal')
     plt.title('Final flux (m$^2$ a$^{-1}$ / 10000)')
@@ -199,7 +202,7 @@ def visualize_eismint2(config, logger, experiment):
     if flwa[timelev, :, :].max() > 0.0:
         # NOT SURE WHICH LEVEL FLWA SHOULD COME FROM - so taking column average
         _contour_mpas(
-            flwa[timelev, :, :].mean(axis=1) * 3600.0*24.0*365.0 / 1.0e-17,
+            flwa[timelev, :, :].mean(axis=1) * 3600.0 * 24.0 * 365.0 / 1.0e-17,
             nCells, xCell, yCell)
     ax.set_aspect('equal')
     # Note: the paper's figure claims units of 10$^{-25}$ Pa$^{-3}$ a$^{-1}$
@@ -215,7 +218,8 @@ def visualize_eismint2(config, logger, experiment):
     # DIVIDE EVOLUTION TIME SERIES
     # ================
     fig = plt.figure(3, facecolor='w')
-    fig.suptitle('Payne et al. Fig. 5, 7, or 8', fontsize=10, fontweight='bold')
+    fig.suptitle('Payne et al. Fig. 5, 7, or 8', fontsize=10,
+                 fontweight='bold')
 
     # get indices for given time
     if experiment == 'b':
@@ -232,14 +236,15 @@ def visualize_eismint2(config, logger, experiment):
     # panel a - thickness
     fig.add_subplot(211)
     timeInd = np.nonzero(years <= endTime)[0][0:]
-    plt.plot(years[timeInd]/1000.0, thickness[timeInd, divideIndex], 'k.-')
+    plt.plot(years[timeInd] / 1000.0, thickness[timeInd, divideIndex], 'k.-')
     plt.ylabel('Thickness (m)')
 
     # panel b - basal temperature
     fig.add_subplot(212)
     # skip the first index cause basalTemperature isn't calculated then
     timeInd = np.nonzero(years <= endTime)[0][1:]
-    plt.plot(years[timeInd]/1000.0, basalTemperature[timeInd, divideIndex], 'k.-')
+    plt.plot(years[timeInd] / 1000.0, basalTemperature[timeInd, divideIndex],
+             'k.-')
     plt.ylabel('Basal temperature (K)')
     plt.xlabel('Time (kyr)')
 
@@ -296,8 +301,8 @@ def visualize_eismint2(config, logger, experiment):
                  'min/mean/max of community', fontsize=10, fontweight='bold')
 
     fig.add_subplot(151)
-    volume = ((thickness[timelev, iceIndices] * areaCell[iceIndices]).sum()
-              / 1000.0**3 / 10.0**6)
+    volume = ((thickness[timelev, iceIndices] * areaCell[iceIndices]).sum() /
+              1000.0**3 / 10.0**6)
     # benchmark results
     plt.plot(np.zeros((3,)), bench['volume'], 'k*')
     if bench['stattype'] == 'relative':
@@ -473,7 +478,7 @@ def _contour_mpas(field, nCells, xCell, yCell, contour_levs=None):
     if len(contour_levs) == 1:
         im = plt.contour(xi, yi, zi)
     else:
-        im = plt.contour(xi, yi, zi, contour_levs, cmap=plt.cm.jet)
+        im = plt.contour(xi, yi, zi, contour_levs, cmap=plt.get_cmap('jet'))
 
     # to see the raw data on top
     # plt.scatter(xCell, yCell, c=temperature[timelev,:,-1], s=100,

compass/landice/tests/ensemble_generator/plot_ensemble.py

@@ -522,7 +522,7 @@ def filter_run():
 
 if plot_maps:
     # axMaps2.plot(GLX, GLY, '.')
-    axMaps2.hist2d(GLX, GLY, (50, 50), cmap=plt.cm.jet)
+    axMaps2.hist2d(GLX, GLY, (50, 50), cmap=plt.get_cmap('jet'))
     figMaps.savefig('figure_maps.png')
 
 if plot_time_series:

compass/ocean/tests/global_convergence/cosine_bell/analysis.py

@@ -131,21 +131,20 @@ def rmse(self, mesh_name):
         # find time since the beginning of run
         ds = xr.open_dataset(f'{mesh_name}_output.nc')
         for j in range(len(ds.xtime)):
-            tt = str(ds.xtime[j].values)
-            tt.rfind('_')
-            DY = float(tt[10:12]) - 1
-            if DY == pd:
+            time_string = str(np.strings.decode(ds.xtime[j].values))
+            day = float(time_string[8:10]) - 1
+            if day == pd:
                 sliceTime = j
                 break
-        HR = float(tt[13:15])
-        MN = float(tt[16:18])
-        t = 86400.0 * DY + HR * 3600. + MN
+        hour = float(time_string[11:13])
+        min = float(time_string[14:16])
+        t = 86400.0 * day + hour * 3600. + min
         # find new location of blob center
         # center is based on equatorial velocity
-        R = init.sphere_radius
-        distTrav = 2.0 * 3.14159265 * R / (86400.0 * pd) * t
+        r = init.sphere_radius
+        distTrav = 2.0 * np.pi * r / (86400.0 * pd) * t
         # distance in radians is
-        distRad = distTrav / R
+        distRad = distTrav / r
         newLon = lonCent + distRad
         if newLon > 2.0 * np.pi:
             newLon -= 2.0 * np.pi
@@ -154,17 +153,17 @@ def rmse(self, mesh_name):
         tracer = np.zeros_like(init.tracer1[0, :, 0].values)
         latC = init.latCell.values
         lonC = init.lonCell.values
-        temp = R * np.arccos(np.sin(latCent) * np.sin(latC) +
-                             np.cos(latCent) * np.cos(latC) * np.cos(
-            lonC - newLon))
+        temp = r * np.arccos(np.sin(latCent) * np.sin(latC) +
+                             (np.cos(latCent) * np.cos(latC) *
+                              np.cos(lonC - newLon)))
         mask = temp < radius
         tracer[mask] = psi0 / 2.0 * (
-            1.0 + np.cos(3.1415926 * temp[mask] / radius))
+            1.0 + np.cos(np.pi * temp[mask] / radius))
 
         # oad forward mode data
         tracerF = ds.tracer1[sliceTime, :, 0].values
         rmseValue = np.sqrt(np.mean((tracerF - tracer)**2))
 
         init.close()
         ds.close()
-        return rmseValue, init.dims['nCells']
+        return rmseValue, init.sizes['nCells']

compass/ocean/tests/sphere_transport/process_output.py

@@ -1,9 +1,10 @@
-import matplotlib.pyplot as plt
-from matplotlib.gridspec import GridSpec
-from matplotlib.colors import ListedColormap
+from importlib import resources
+
 import matplotlib
+import matplotlib.pyplot as plt
 import numpy as np
-from importlib import resources
+from matplotlib.colors import ListedColormap
+from matplotlib.gridspec import GridSpec
 
 
 def appx_mesh_size(dataset):
@@ -260,8 +261,8 @@ def read_ncl_rgb_file(cmap_filename):
     map_file_found = False
     try:
         with resources.open_text(
-                "compass.ocean.tests.sphere_transport.resources", cmap_filename) \
-                as f:
+                "compass.ocean.tests.sphere_transport.resources",
+                cmap_filename) as f:
             flines = f.readlines()
         map_file_found = True
     except BaseException:
@@ -277,7 +278,7 @@ def read_ncl_rgb_file(cmap_filename):
         result = ListedColormap(rgb, name=cmap_filename)
     else:
         print("error reading ncl colormap. using matplotlib default instead.")
-        result = matplotlib.cm.get_cmap()
+        result = matplotlib.pyplot.get_cmap()
     return result
 
 
@@ -364,8 +365,7 @@ def plot_sol(fig, tcname, dataset):
                              cmap="seismic",
                              vmin=diffmin,
                              vmax=diffmax)
-    lcm = axes[9].tricontourf(xc, yc, dataset.variables["layerThickness"]
-                              [0, :, 1])
+    axes[9].tricontourf(xc, yc, dataset.variables["layerThickness"][0, :, 1])
     axes[9].set_ylabel('layer thickness')
     axes[10].tricontourf(xc, yc, dataset.variables["layerThickness"]
                          [0, :, 1])
@@ -389,8 +389,8 @@ def plot_sol(fig, tcname, dataset):
         axes[9 + i].set_xticklabels(xticklabels)
     for i in range(9):
         axes[i].set_xticklabels([])
-    cb1 = fig.colorbar(cm, ax=axes[8])
-    cb2 = fig.colorbar(tcm, ax=axes[5])
+    fig.colorbar(cm, ax=axes[8])
+    fig.colorbar(tcm, ax=axes[5])
     # cb3 = fig.colorbar(lcm, ax=axes[11])
     fig.suptitle(tcname)