Add a tides test case for the VR45to5 mesh

compass/ocean/tests/tides/__init__.py

@@ -1,8 +1,7 @@
-from compass.testgroup import TestGroup
-
-from compass.ocean.tests.tides.mesh import Mesh
-from compass.ocean.tests.tides.init import Init
 from compass.ocean.tests.tides.forward import Forward
+from compass.ocean.tests.tides.init import Init
+from compass.ocean.tests.tides.mesh import Mesh
+from compass.testgroup import TestGroup
 
 
 class Tides(TestGroup):
@@ -17,7 +16,7 @@ def __init__(self, mpas_core):
         super().__init__(mpas_core=mpas_core,
                          name='tides')
 
-        for mesh_name in ['Icos7']:
+        for mesh_name in ['Icos7', 'VR45to5']:
 
             mesh = Mesh(test_group=self, mesh_name=mesh_name)
             self.add_test_case(mesh)

compass/ocean/tests/tides/analysis/__init__.py

@@ -1,14 +1,15 @@
-from compass.step import Step
+import os
 
-import netCDF4
-import matplotlib.pyplot as plt
-import matplotlib.cm as cm
-import numpy as np
 import cartopy.crs as ccrs
 import cartopy.feature as cfeature
-import os
+import matplotlib.cm as cm
+import matplotlib.pyplot as plt
+import netCDF4
+import numpy as np
 from mpas_tools.logging import check_call
 
+from compass.step import Step
+
 
 class Analysis(Step):
     """
@@ -103,7 +104,7 @@ def write_coordinate_file(self, idx):
         # Write coordinate file for OTPS2
         f = open('lat_lon', 'w')
         for i in range(nCells):
-            f.write(str(lat_grid[i])+'  '+str(lon_grid[i])+'\n')
+            f.write(str(lat_grid[i]) + '  ' + str(lon_grid[i]) + '\n')
         f.close()
 
     def setup_otps2(self):
@@ -121,7 +122,7 @@ def setup_otps2(self):
                       'comment': '! 2. latitude/longitude/<time> file'},
                      {'inp': 'z',
                       'comment': '! 3. z/U/V/u/v'},
-                     {'inp':  con,
+                     {'inp': con,
                       'comment': '! 4. tidal constituents to include'},
                      {'inp': 'AP',
                       'comment': '! 5. AP/RI'},
@@ -137,10 +138,10 @@ def setup_otps2(self):
                 os.mkdir('inputs')
 
             # Write the setup_con file
-            f = open('inputs/'+con+'_setup', 'w')
+            f = open(f'inputs/{con}_setup', 'w')
             for line in lines:
                 spaces = 28 - len(line['inp'])
-                f.write(line['inp'] + spaces*' ' + line['comment'] + '\n')
+                f.write(line['inp'] + spaces * ' ' + line['comment'] + '\n')
             f.close()
 
             # Write the Model_atlas_con file
@@ -180,18 +181,18 @@ def read_otps2_output(self, idx):
                 line_sp = line.split()
                 if line_sp[2] != '*************':
                     val = float(line_sp[2])
-                    self.mesh_AP[con]['amp'][start+i] = val
+                    self.mesh_AP[con]['amp'][start + i] = val
                 else:
-                    self.mesh_AP[con]['amp'][start+i] = -9999
+                    self.mesh_AP[con]['amp'][start + i] = -9999
 
                 if line_sp[3] != 'Site':
                     val = float(line_sp[3])
                     if val < 0:
                         val = val + 360.0
-                    self.mesh_AP[con]['phase'][start+i] = val
+                    self.mesh_AP[con]['phase'][start + i] = val
 
                 else:
-                    self.mesh_AP[con]['phase'][start+i] = -9999
+                    self.mesh_AP[con]['phase'][start + i] = -9999
 
     def append_tpxo_data(self):
         """
@@ -268,8 +269,8 @@ def plot(self):
 
         lon_grid = np.mod(data_nc.variables['lonCell'][:] + np.pi,
                           2.0 * np.pi) - np.pi
-        lon_grid = lon_grid*180.0/np.pi
-        lat_grid = data_nc.variables['latCell'][:]*180.0/np.pi
+        lon_grid = lon_grid * 180.0 / np.pi
+        lat_grid = data_nc.variables['latCell'][:] * 180.0 / np.pi
 
         nCells = lon_grid.size
         data1 = np.zeros((nCells))
@@ -287,14 +288,14 @@ def plot(self):
         # Use these to fix up the plots
         subplot_ticks = [[np.linspace(0, 0.65, 10), np.linspace(0, 0.65, 10),
                           np.linspace(0, 0.13, 10), np.linspace(0, 0.13, 10)],
-                         [np.linspace(0, 1.4,  10), np.linspace(0, 1.4,  10),
+                         [np.linspace(0, 1.40, 10), np.linspace(0, 1.40, 10),
                           np.linspace(0, 0.22, 10), np.linspace(0, 0.25, 10)],
                          [np.linspace(0, 0.22, 10), np.linspace(0, 0.22, 10),
                           np.linspace(0, 0.05, 10), np.linspace(0, 0.05, 10)],
-                         [np.linspace(0, 0.5,  10), np.linspace(0, 0.5,  10),
+                         [np.linspace(0, 0.50, 10), np.linspace(0, 0.50, 10),
                           np.linspace(0, 0.08, 10), np.linspace(0, 0.08, 10)],
-                         [np.linspace(0, 0.7,  10), np.linspace(0, 0.7,  10),
-                          np.linspace(0, 0.5,  10), np.linspace(0, 0.5,  10)]]
+                         [np.linspace(0, 0.70, 10), np.linspace(0, 0.70, 10),
+                          np.linspace(0, 0.50, 10), np.linspace(0, 0.50, 10)]]
 
         for i, con in enumerate(constituent_list):
 
@@ -311,39 +312,45 @@ def plot(self):
             data2_phase[:] = data_nc.variables[
                 f'{con}Phase{self.tpxo_version}'][:]
 
-            data1_phase = data1_phase*np.pi/180.0
-            data2_phase = data2_phase*np.pi/180.0
+            data1_phase = data1_phase * np.pi / 180.0
+            data2_phase = data2_phase * np.pi / 180.0
 
             # Calculate RMSE values
-            rmse_amp = 0.5*(data1 - data2)**2
-            rmse_com = 0.5*(data2**2 + data1**2) \
-                - data1*data2*np.cos(data2_phase - data1_phase)
+            rmse_amp = 0.5 * (data1 - data2)**2
+            rmse_com = 0.5 * (data2**2 + data1**2) \
+                - data1 * data2 * np.cos(data2_phase - data1_phase)
 
             # Calculate mean (global) values
-            idx = np.where((depth > 20)
-                           & (rmse_com < 1000) & (rmse_amp < 1000))
+            idx = np.where((depth > 20) &
+                           (rmse_com < 1000) & (rmse_amp < 1000))
             area_tot = np.sum(area[idx])
-            glo_rmse_amp = np.sqrt(np.sum(rmse_amp[idx]*area[idx])/area_tot)
-            glo_rmse_com = np.sqrt(np.sum(rmse_com[idx]*area[idx])/area_tot)
+            glo_rmse_amp = np.sqrt(np.sum(rmse_amp[idx] * area[idx]) /
+                                   area_tot)
+            glo_rmse_com = np.sqrt(np.sum(rmse_com[idx] * area[idx]) /
+                                   area_tot)
             print('Global RMSE (Amp) = ', glo_rmse_amp)
             print('Global RMSE (Com) = ', glo_rmse_com)
 
             # Calculate shallow RMSE (<=1000m)
-            idx = np.where((depth > 20) & (depth < 1000)
-                           & (np.abs(lat_grid) < 66)
-                           & (rmse_com < 1000) & (rmse_amp < 1000))
+            idx = np.where((depth > 20) & (depth < 1000) &
+                           (np.abs(lat_grid) < 66) &
+                           (rmse_com < 1000) & (rmse_amp < 1000))
             area_tot = np.sum(area[idx])
-            shal_rmse_amp = np.sqrt(np.sum(rmse_amp[idx]*area[idx])/area_tot)
-            shal_rmse_com = np.sqrt(np.sum(rmse_com[idx]*area[idx])/area_tot)
+            shal_rmse_amp = np.sqrt(np.sum(rmse_amp[idx] * area[idx]) /
+                                    area_tot)
+            shal_rmse_com = np.sqrt(np.sum(rmse_com[idx] * area[idx]) /
+                                    area_tot)
             print('Shallow RMSE (Amp) = ', shal_rmse_amp)
             print('Shallow RMSE (Com) = ', shal_rmse_com)
 
             # Calculate deep RMSE (>1000m)
-            idx = np.where((depth >= 1000) & (np.abs(lat_grid) < 66)
-                           & (rmse_com < 1000) & (rmse_amp < 1000))
+            idx = np.where((depth >= 1000) & (np.abs(lat_grid) < 66) &
+                           (rmse_com < 1000) & (rmse_amp < 1000))
             area_tot = np.sum(area[idx])
-            deep_rmse_amp = np.sqrt(np.sum(rmse_amp[idx]*area[idx])/area_tot)
-            deep_rmse_com = np.sqrt(np.sum(rmse_com[idx]*area[idx])/area_tot)
+            deep_rmse_amp = np.sqrt(np.sum(rmse_amp[idx] * area[idx]) /
+                                    area_tot)
+            deep_rmse_com = np.sqrt(np.sum(rmse_com[idx] * area[idx]) /
+                                    area_tot)
             print('Deep RMSE (Amp) = ', deep_rmse_amp)
             print('Deep RMSE (Com) = ', deep_rmse_com)
 
@@ -353,12 +360,12 @@ def plot(self):
             # Plot data
             fig = plt.figure(figsize=(18, 12))
             subplot_title = [f'{con} Amplitude (simulation) [m]',
-                             f'{con} Amplitude (TPXO8) [m]',
+                             f'{con} Amplitude ({self.tpxo_version}) [m]',
                              f'{con} RMSE (Amplitude) [m]',
                              f'{con} RMSE (Complex) [m]']
 
             for subplot in range(0, 4):
-                ax = fig.add_subplot(2, 2, subplot+1,
+                ax = fig.add_subplot(2, 2, subplot + 1,
                                      projection=ccrs.PlateCarree())
                 ax.set_title(subplot_title[subplot], fontsize=20)
                 levels = subplot_ticks[i][subplot][:]
@@ -407,9 +414,9 @@ def plot(self):
                 cbar.ax.tick_params(labelsize=16)
 
             fig.tight_layout()
-            global_err = str(round(glo_rmse_com*100, 3))
-            deep_err = str(round(deep_rmse_com*100, 3))
-            shallow_err = str(round(shal_rmse_com*100, 3))
+            global_err = str(round(glo_rmse_com * 100, 3))
+            deep_err = str(round(deep_rmse_com * 100, 3))
+            shallow_err = str(round(shal_rmse_com * 100, 3))
             fig.suptitle(f'Complex RMSE: Global = {global_err} cm; '
                          f'Deep = {deep_err} cm; '
                          f'Shallow = {shallow_err} cm',
@@ -430,7 +437,7 @@ def run(self):
 
         # Setup chunking for TPXO extraction with large meshes
         indices = np.arange(self.nCells)
-        nchunks = np.ceil(self.nCells/200000)
+        nchunks = np.ceil(self.nCells / 200000)
         index_chunks = np.array_split(indices, nchunks)
 
         # Initialize data structure for TPXO values

compass/ocean/tests/tides/dem/__init__.py

@@ -0,0 +1,59 @@
+import os
+
+import compass.ocean.tests.tides.dem.dem_pixel as dem_pixel
+from compass.step import Step
+
+
+class CreatePixelFile(Step):
+    """
+    A step for creating a pixel file for creating MPAS meshes
+
+    Attributes
+    ----------
+    """
+
+    def __init__(self, test_case):
+        """
+        Create the step
+
+        Parameters
+        ----------
+        test_case : compass.ocean.tests.tides.init.Init
+            The test case this step belongs to
+        """
+        super().__init__(test_case=test_case, name='pixel',
+                         ntasks=1, min_tasks=1, openmp_threads=1)
+
+        self.add_input_file(
+            filename='GEBCO_2023_sub_ice_topo.nc',
+            target='GEBCO_2023_sub_ice_topo.nc',
+            database='bathymetry_database')
+
+        self.add_input_file(
+            filename='RTopo-2.0.4_30sec_bedrock_topography.nc',
+            target='RTopo-2.0.4_30sec_bedrock_topography.nc',
+            database='bathymetry_database')
+
+        self.add_input_file(
+            filename='RTopo-2.0.4_30sec_surface_elevation.nc',
+            target='RTopo-2.0.4_30sec_surface_elevation.nc',
+            database='bathymetry_database')
+
+        self.add_input_file(
+            filename='RTopo-2.0.4_30sec_ice_base_topography.nc',
+            target='RTopo-2.0.4_30sec_ice_base_topography.nc',
+            database='bathymetry_database')
+
+    def run(self):
+        """
+        Run this step of the test case
+        """
+
+        self.init_path = './'
+
+        if not os.path.exists('RTopo_2_0_4_30sec_pixel.nc'):
+            dem_pixel.rtopo_30sec(self.init_path, self.init_path)
+        if not os.path.exists('GEBCO_v2023_30sec_pixel.nc'):
+            dem_pixel.gebco_30sec(self.init_path, self.init_path)
+        if not os.path.exists('RTopo_2_0_4_GEBCO_v2023_30sec_pixel.nc'):
+            dem_pixel.rtopo_gebco_30sec(self.init_path, self.init_path)

compass/ocean/tests/tides/dem/dem_names.py

@@ -0,0 +1,21 @@
+
+# -- long_name strings for NetCDF variables
+
+class base:
+    pass
+
+
+names = base()
+names.bed_elevation = "elevation of bed"
+names.ocn_thickness = "thickness of ocn"
+names.ice_thickness = "thickness of ice"
+names.ocn_cover = "fractional cover of ocn, 0-1"
+names.ice_cover = "fractional cover of ice, 0-1"
+names.bed_slope = "RMS magnitude of bed slopes"
+names.bed_slope_deg = "arc-tangent of RMS bed slopes"
+names.bed_dz_dx = "derivative of bed elevation along lon.-axis"
+names.bed_dz_dy = "derivative of bed elevation along lat.-axis"
+names.bed_elevation_profile = "sub-grid percentiles of bed elevation"
+names.bed_slope_profile = "sub-grid percentiles of RMS bed slope"
+names.ocn_thickness_profile = "sub-grid percentiles of ocn thickness"
+names.ice_thickness_profile = "sub-grid percentiles of ice thickness"