dev

davidhassell · davidhassell · commit bc279f2bdb26 · 2023-08-23T18:31:14.000+01:00
diff --git a/cf/field.py b/cf/field.py
@@ -2365,55 +2365,41 @@ def _weights_geometry_area(
 
             print (111)
             
-#            if ugrid:
-#                # UGRID face cell bounds do not have a "parts"
-#                # dimension, so insert a dummy one.
-#                 x = x.insert_dimension(1).persist()
-#                y = y.insert_dimension(1).persist()
-#                x = x.concatenate(x, x[:,0:1]], axis=1)
-#                y = y.concatenate(y, y[:,0:1]], axis=1)
+            if ugrid:
+                x = x.array
+                col_indices = np.ma.count(x, axis=1)
+                empty_column = np.ma.masked_all((x.shape[0], 1), dtype=x.dtype)
+                x = np.ma.concatenate((x, empty_column), axis=1)
+                x[np.arange(x.shape[0]), col_indices] = x[:, 0]
+#                x = np.ma.where(x <0, x + np.pi, x)
+                x = self._Data(x, units= _units_radians)
+                print (x.shape, x[0].array)
+                
+                y = y.array
+                col_indices = np.ma.count(y, axis=1)
+                empty_column = np.ma.masked_all((y.shape[0], 1), dtype=y.dtype)
+                y = np.ma.concatenate((y, empty_column), axis=1)
+                y[np.arange(y.shape[0]), col_indices] = y[:, 0]
+                y = self._Data(y, units= _units_radians)
+                print(y.shape, y[0].array)
 
-
-            
             interior_angle = self._weights_interior_angle(x, y)
-
+            print (interior_angle.array[0])
             
             # Find the number of edges of each polygon (note that
             # this number may be one too few, but we'll adjust for
             # that later).
 #            N = interior_angle.sample_size(-1, squeeze=True)
             N = interior_angle.count(-1, keepdims=False)
+            print (N.array)
 
             all_areas = interior_angle.sum(-1, squeeze=True) - (N - 2) * np.pi
-
-            if ugrid:                
-                x = x.persist()
-                y = y.persist()
-                for i, (nodes_x, nodes_y) in enumerate(zip(x, y)):
-                    print (axis, i)
-                    nodes_x = nodes_x.compressed().persist()
-                    nodes_y = nodes_y.compressed().persist()
-                    
-                    if (nodes_x.size and nodes_x[0] != nodes_x[-1]) or (
-                            nodes_y.size and nodes_y[0] != nodes_y[-1]
-                    ):
-                        # First and last nodes of this polygon
-                        # part are different => need to account
-                        # for the "last" edge of the polygon that
-                        # joins the first and last points.
-                        interior_angle = self._weights_interior_angle(
-                            nodes_x[[0, -1]], nodes_y[[0, -1]]
-                        )
-                        all_areas[i] += interior_angle + np.pi
-
-            else:
+            print (all_areas.array[0])
+            if not ugrid:                
                 for i, (parts_x, parts_y) in enumerate(zip(x, y)):
-                    print (ugrid, axis, i, parts_x.shape, parts_y.shape)
                     for j, (nodes_x, nodes_y) in enumerate(zip(parts_x, parts_y)):
                         nodes_x = nodes_x.compressed()
                         nodes_y = nodes_y.compressed()
-                        print (nodes_x.shape)
-    
                         if (nodes_x.size and nodes_x[0] != nodes_x[-1]) or (
                             nodes_y.size and nodes_y[0] != nodes_y[-1]
                         ):
@@ -2425,7 +2411,7 @@ def _weights_geometry_area(
                                 nodes_x[[0, -1]], nodes_y[[0, -1]]
                             )
     
-                            all_areas[i, j] += interior_angle + np.pi
+                            all_areas[i, j] += interior_angle + np.pi # -?
 
             area_min = all_areas.min()
             if area_min < 0:
@@ -2446,7 +2432,9 @@ def _weights_geometry_area(
         if ugrid:
             areas = all_areas
         else:
-            areas = all_areas.sum(-1, squeeze=True)
+            # Sum the areas of each geometry polygon part to get the
+            # total area of each cell
+             areas = all_areas.sum(-1, squeeze=True)
 
         if measure and spherical and aux_Z is not None:
             # Multiply by radius squared, accounting for any Z
@@ -2787,7 +2775,7 @@ def _weights_interior_angle(self, data_lambda, data_phi):
                 The interior angles in units of radians.
 
         """
-        delta_lambda = data_lambda.diff(axis=-1)
+        delta_lambda = abs(data_lambda.diff(axis=-1))
 
         cos_phi = data_phi.cos()
         sin_phi = data_phi.sin()
@@ -2813,6 +2801,7 @@ def _weights_interior_angle(self, data_lambda, data_phi):
         denominator = (
             sin_phi_1 * sin_phi_2 + cos_phi_1 * cos_phi_2 * cos_delta_lambda
         )
+        print ('dem=', denominator[0].array)
 
         # TODO RuntimeWarning: overflow encountered in true_divide comes from
         # numerator/denominator with missing values
