add new NEMO reader

PyFVCOM/grid/_grid.py

@@ -1703,6 +1703,12 @@ class OpenBoundary(object):
             # Keep positive down depths.
             z = -self.sigma.layers_z
 
+        ## 2D lat lon to 1D lat, lon  , maz (is this okay?)
+        #if len(coarse.grid.lat)==2:
+        #    coarse.grid.lat = coarse.grid.lat[:,0]
+        #if len(coarse.grid.lon)==2:
+        #    coarse.grid.lon = coarse.grid.lon[0,:]
+
         if constrain_coordinates:
             x[x < coarse.grid.lon.min()] = coarse.grid.lon.min()
             x[x > coarse.grid.lon.max()] = coarse.grid.lon.max()
@@ -1713,13 +1719,16 @@ class OpenBoundary(object):
             # move it to the nearest in grid point if so.
             if not mode == 'surface':
                 land_mask = getattr(coarse.data, coarse_name)[0, ...].mask[0, :, :]
+                #land_mask = coarse.tmask[0, 0, ...]     # maz
             else:
                 land_mask = getattr(coarse.data, coarse_name)[0, ...].mask
+                #land_mask = coarse.tmask[0, ...]        # maz
 
             sea_points = np.ones(land_mask.shape)
             sea_points[land_mask] = np.nan
 
             ft_sea = RegularGridInterpolator((coarse.grid.lat, coarse.grid.lon), sea_points, method='linear', fill_value=np.nan)
+#            ft_sea = RegularGridInterpolator((coarse.grid.nav_lat[0,:], coarse.grid.nav_lon[:,0]), sea_points, method='linear', fill_value=np.nan)
             internal_points = np.isnan(ft_sea(np.asarray([y, x]).T))
 
             if np.any(internal_points):

PyFVCOM/preproc.py

@@ -5019,7 +5019,7 @@ class RegularReader(FileReader):
             return np.unravel_index(index, self.grid.lon.shape)
 
 
-class NEMOReader(RegularReader):
+class NEMOReaderNew(RegularReader):
     """
     Class to read in regularly gridded model output. This provides a similar interface to a PyFVCOM.read.FileReader
     object but with an extra spatial dimension. This is currently based on CMEMS model outputs (i.e. NEMO).
@@ -5027,7 +5027,7 @@ class NEMOReader(RegularReader):
     Author(s)
     ---------
     Pierre Cazenave (Plymouth Marine Laboratory)
-
+    
     Credits
     -------
     This code leverages ideas (and in some cases, code) from PySeidon (https://github.com/GrumpyNounours/PySeidon)
@@ -5082,12 +5082,16 @@ class NEMOReader(RegularReader):
                 if 'deptht' in self._dims:
                     self.tmask = self.tmask[..., self._dims['deptht'], :, :]
                 # Make it so it fits with time-varying data.
-                self.tmask = np.tile(np.squeeze(self.tmask), [self.dims.time_counter, 1, 1, 1])
+                #self.tmask = np.tile(np.squeeze(self.tmask), [self.dims.time_counter, 1, 1, 1])
+                self.tmask = np.tile(np.squeeze(self.tmask), [self.dims.time, 1, 1, 1])
             for var in self.data:
                 # We could use masking here, but this feels more bulletproof (I think some bits of numpy/scipy ignore
                 # masks when interpolating).
                 current_data = getattr(self.data, var)
-                current_data[~self.tmask] = np.nan
+                #current_data[~self.tmask] = np.nan
+                if len(current_data.shape) == 3:
+                    self.tmask = self.tmask[:,0,:,:]
+                current_data[self.tmask] = np.nan        
                 setattr(self.data, var, current_data)
         else:
             # If we don't have a tmask file, we'll try our best to minimise potential issues with crappy NEMO data.
@@ -5117,7 +5121,7 @@ class NEMOReader(RegularReader):
         if hasattr(self.dims, 'lat'):
             yname = 'lat'
             ydim = self.dims.lat
-        elif hasattr(self.dims, 'x'):
+        elif hasattr(self.dims, 'y'):
             yname = 'y'
             ydim = self.dims.y
         else:
@@ -5198,8 +5202,11 @@ class NEMOReader(RegularReader):
 
         """
         if grid_variables is None:
-            grid_variables = {'lon': 'nav_lon', 'nav_lat': 'lat', 'x': 'x', 'y': 'y', 'depth': 'depth',
+            grid_variables = {'lon': 'nav_lon', 'lat': 'nav_lat', 'x': 'x', 'y': 'y', 'depth': 'deptht',
                               'Longitude': 'Longitude', 'Latitude': 'Latitude'}
+            self.dims.lon = self.dims.x    
+            self.dims.lat = self.dims.y   
+
 
         self.grid = PassiveStore()
         # Get the grid data.
@@ -5265,6 +5272,7 @@ class NEMOReader(RegularReader):
         self.grid.lon = np.unique(self.grid.nav_lon)
         self.grid.lat = np.unique(self.grid.nav_lat)
 
+
         # Check if we've been given vertical dimensions to subset in too, and if so, do that. Check we haven't
         # already done this if the 'node' and 'nele' sections above first.
         for var in ['depth']: