diff --git a/conda_package/mpas_tools/ocean/transects.py b/conda_package/mpas_tools/ocean/transects.py index a8940a971..bf6a140f3 100644 --- a/conda_package/mpas_tools/ocean/transects.py +++ b/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 diff --git a/conda_package/mpas_tools/viz/transects.py b/conda_package/mpas_tools/viz/transects.py index c592aff1b..0d9670b47 100644 --- a/conda_package/mpas_tools/viz/transects.py +++ b/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