Updates(vtk): updates for errant interpolation in point_scalar routin…

…es (#1601)

* updates to UnstructuredGrid and VertexGrid verts property. Rotate and translate verts based on modelgrid offsets. These has been updated for consistency in behavior with StructuredGrid

flopy/discretization/unstructuredgrid.py

@@ -5,8 +5,7 @@
 import numpy as np
 from matplotlib.path import Path
 
-from ..utils.geometry import is_clockwise
-from ..utils.gridutil import get_lni
+from ..utils.geometry import is_clockwise, transform
 from .grid import CachedData, Grid
 
 
@@ -225,7 +224,17 @@ def verts(self):
         if self._vertices is None:
             return self._vertices
         else:
-            return np.array([list(t)[1:] for t in self._vertices], dtype=float)
+            verts = np.array(
+                [list(t)[1:] for t in self._vertices], dtype=float
+            ).T
+            x, y = transform(
+                verts[0],
+                verts[1],
+                self.xoffset,
+                self.yoffset,
+                self.angrot_radians,
+            )
+            return np.array(list(zip(x, y)))
 
     @property
     def iac(self):

flopy/discretization/vertexgrid.py

@@ -5,7 +5,7 @@
 import numpy as np
 from matplotlib.path import Path
 
-from ..utils.geometry import is_clockwise
+from ..utils.geometry import is_clockwise, transform
 from .grid import CachedData, Grid
 
 
@@ -147,7 +147,11 @@ def cell2d(self):
 
     @property
     def verts(self):
-        return np.array([list(t)[1:] for t in self._vertices], dtype=float)
+        verts = np.array([list(t)[1:] for t in self._vertices], dtype=float).T
+        x, y = transform(
+            verts[0], verts[1], self.xoffset, self.yoffset, self.angrot_radians
+        )
+        return np.array(list(zip(x, y)))
 
     @property
     def shape(self):

flopy/export/vtk.py

@@ -353,9 +353,10 @@ def _create_point_scalar_graph(self):
         for k, d in graph.items():
             for _, pt in enumerate(d["vtk_points"]):
                 for ixx, value in enumerate(d["idx"]):
-                    numpy_graph[ixx, pt] = value
-                    xvert[ixx, pt] = d["xv"][ixx]
-                    yvert[ixx, pt] = d["yv"][ixx]
+                    if self.idomain[value] > 0:
+                        numpy_graph[ixx, pt] = value
+                        xvert[ixx, pt] = d["xv"][ixx]
+                        yvert[ixx, pt] = d["yv"][ixx]
 
         # now create the IDW weights for point scalars
         xc = np.ravel(self.modelgrid.xcellcenters)
@@ -647,16 +648,22 @@ def _build_point_scalar_array(self, array):
                 for ix, pt in enumerate(value["vtk_points"]):
                     ps_array[pt] = array[value["idx"][ix]]
         else:
-            ps_graph = self._point_scalar_numpy_graph
+            ps_graph = self._point_scalar_numpy_graph.copy()
+            idxs = np.where(np.isnan(array))
+            not_graphed = np.isin(ps_graph, idxs[0])
+            ps_graph[not_graphed] = -1
             ps_array = np.where(ps_graph >= 0, array[ps_graph], np.nan)
 
             # do inverse distance weighting and apply mask to retain
             # nan valued cells because numpy returns 0 when all vals are nan
-            weighted_vals = self._idw_weight_graph * ps_array
+            weight_graph = self._idw_weight_graph.copy()
+            weight_graph[not_graphed] = np.nan
+            weighted_vals = weight_graph * ps_array
             mask = np.isnan(weighted_vals).all(axis=0)
             weighted_vals = np.nansum(weighted_vals, axis=0)
             weighted_vals[mask] = np.nan
-            ps_array = weighted_vals / self._idw_total_weight_graph
+            total_weight_graph = np.nansum(weight_graph, axis=0)
+            ps_array = weighted_vals / total_weight_graph
 
         return ps_array