Merge pull request #534 from xylar/fix-planar-transect

Fixes to transect extraction and plotting

conda_package/mpas_tools/ocean/transects.py

@@ -378,9 +378,14 @@ def _get_vertical_coordinate(dsTransect, layerThickness, bottomDepth,
     layerThicknessTransect = (layerThicknessTransect*interpCellWeights).sum(
         dim='nHorizWeights')
 
+    interpMask = maxLevelCell.isel(nCells=interpCellIndices) >= 0
+    interpHorizCellWeights = interpMask*dsTransect.interpHorizCellWeights
+    weightSum = interpHorizCellWeights.sum(dim='nHorizWeights')
+    interpHorizCellWeights = \
+        (interpHorizCellWeights/weightSum).where(interpMask)
+
     sshTransect = ssh.isel(nCells=interpCellIndices)
-    sshTransect = (sshTransect*dsTransect.interpHorizCellWeights).sum(
-        dim='nHorizWeights')
+    sshTransect = (sshTransect*interpHorizCellWeights).sum(dim='nHorizWeights')
 
     zBot = sshTransect - layerThicknessTransect.cumsum(dim='nVertLevels')
     zTop = zBot + layerThicknessTransect

conda_package/mpas_tools/viz/transects.py

@@ -1,6 +1,6 @@
+import matplotlib.pyplot as plt
 import numpy
 import xarray
-
 from scipy.spatial import cKDTree
 from shapely.geometry import LineString, Point
 
@@ -516,9 +516,9 @@ def find_planar_transect_cells_and_weights(xTransect, yTransect, dsTris, dsMesh,
         xIntersection = numpy.array(xIntersection)
         yIntersection = numpy.array(yIntersection)
         nodeWeights = numpy.array(nodeWeights)
-        node0Inter = numpy.array(node0Inter)
-        node1Inter = numpy.array(node1Inter)
-        trisInter = numpy.array(trisInter)
+        node0Inter = numpy.array(node0Inter, dtype=int)
+        node1Inter = numpy.array(node1Inter, dtype=int)
+        trisInter = numpy.array(trisInter, dtype=int)
 
         dNodeLocal = dStart + distances
 
@@ -626,6 +626,17 @@ def _sort_intersections(dNode, tris, nodes, xOut, yOut, zOut, interpCells,
     dSorted = dNode[sortIndices]
     trisSorted = tris[sortIndices]
 
+    # sometimes we end up with redundant intersections in the transect, and
+    # these need to be removed
+    unique_d_tris = dict()
+    for index in range(len(dSorted)):
+        unique_d_tris[(dSorted[index], trisSorted[index])] = index
+
+    unique_indices = list(unique_d_tris.values())
+    sortIndices = sortIndices[unique_indices]
+    dSorted = dSorted[unique_indices]
+    trisSorted = trisSorted[unique_indices]
+
     nodesAreSame = numpy.abs(dSorted[1:] - dSorted[:-1]) < epsilon
     if nodesAreSame[0]:
         # the first two nodes are the same, so the first transect point is in a