dev

cf/data/data.py

@@ -178,6 +178,7 @@ def __init__(
         copy=True,
         dtype=None,
         mask=None,
+            masked_values=None,
         to_memory=False,
         init_options=None,
         _use_array=True,
@@ -484,6 +485,10 @@ def __init__(
         if mask is not None:
             self.where(mask, cf_masked, inplace=True)
 
+        # Apply masked values
+        if masked_values is not None:
+            self.masked_values(masked_values, inplace=True)
+
     @property
     def dask_compressed_array(self):
         """Returns a dask array of the compressed data.

cf/weights.py

@@ -574,8 +574,9 @@ def polygon_area(
         else:
             y.Units = x.Units
 
-        cols = x.count(axis=-1).array
-        if (y.count(axis=-1) != cols).any():
+        cols = x.count(axis=-1, keepdims=False).array
+        print(repr(cols), cols.shape) #, cols[:, 0], cols[:,1])
+        if (y.count(axis=-1, keepdims=False) != cols).any():
             raise ValueError(
                 "Can't create area weights for "
                 f"{f.constructs.domain_axis_identity(axis)!r} axis: "
@@ -598,10 +599,22 @@ def polygon_area(
             # Ascertain whether or not the first and last boundary
             # vertices are coincident in all polygons. Note that this
             # is never the case for UGRID cells.
+            duplicate = np.array(False)
+        else:
             cols_m1 = cols - 1
-            duplicate = x[:, ..., 0].isclose(x[rows, ..., cols_m1]) & y[
-                :, ..., 0
-            ].isclose(y[rows, ..., cols_m1])
+            if x.ndim == 3:
+                ir = 0
+            else:
+                ir= Ellipsis
+
+            if equal_nodes:
+                c = cols_m1
+            else:
+                c = cols_m1[:, 0]
+
+            duplicate = x[:, ir, 0].isclose(x[rows, ir, cols_m1]) & y[
+                :, ir, 0
+            ].isclose(y[rows, ir, cols_m1])
             duplicate = duplicate.array
 
         x = x.array
@@ -611,7 +624,10 @@ def polygon_area(
         # and last bounds neighbours
         empty_col_x = np.ma.masked_all(x.shape[:-1] + (1,), dtype=x.dtype)
         empty_col_y = np.ma.masked_all(y.shape[:-1] + (1,), dtype=y.dtype)
+        print ('duplicate=', duplicate)
 
+        print ('x=',x)
+        print ('y=',y)
         if ugrid or not duplicate.any():
             # The first and last boundary vertices are different in
             # all polygons, i.e. No. nodes = No. edges.
@@ -624,13 +640,17 @@ def polygon_area(
             y = np.ma.concatenate((empty_col_y, y, empty_col_y), axis=-1)
 
             if x.ndim == 3 and not equal_nodes:
+                print (99999)
                 for i in range(x.shape[1]):
                     c = cols[:, i]
                     x[:, i, 0] = x[rows, i, c]
                     y[:, i, 0] = y[rows, i, c]
                     c = cols_p1[:, i]
                     x[rows, i, c] = x[rows, i, 1]
                     y[rows, i, c] = y[rows, i, 1]
+
+                print ('x=',x)
+                print ('y=',y)
             else:
                 x[:, ..., 0] = x[rows, ..., cols]
                 y[:, ..., 0] = y[rows, ..., cols]
@@ -681,6 +701,7 @@ def polygon_area(
             # (https://en.wikipedia.org/wiki/Spherical_trigonometry).
             # --------------------------------------------------------
             interior_angles = cls.interior_angles(f, x, y, interior_ring)
+            print ('interior_angles=',interior_angles.array)
             all_areas = interior_angles.sum(-1, squeeze=True) - (N - 2) * np.pi
         else:
             # --------------------------------------------------------
@@ -750,6 +771,8 @@ def polygon_area(
         if return_areas:
             return areas
 
+        print ('areas=', areas)
+
         weights[(axis,)] = areas
         weights_axes.add(axis)