import argparse
import logging
import sys
import numpy
import xarray
from geometric_features.aggregation.ocean import moc
import mpas_tools.mesh.conversion
from mpas_tools.io import default_nchar, write_netcdf
[docs]
def make_moc_basins_and_transects(
gf,
mesh_filename,
mask_and_transect_filename,
geojson_filename=None,
mask_filename=None,
logger=None,
dir=None,
):
"""
Builds features defining the ocean basins and southern transects used in
computing the meridional overturning circulation (MOC)
Parameters
----------
gf : geometric_features.GeometricFeatures
An object that knows how to download and read geometric features
mesh_filename : str
A file with MPAS mesh information
mask_and_transect_filename : str
A file to write the MOC region masks and southern-boundary transects to
geojson_filename : str, optional
A file to write MOC regions to
mask_filename : str, optional
A file to write MOC region masks to
logger : ``logging.Logger``, optional
A logger for the output if not stdout
dir : str, optional
A directory in which a temporary directory will be added with files
produced during conversion and then deleted upon completion.
Returns
-------
fc : geometric_features.FeatureCollection
The new feature collection
"""
# Authors
# -------
# Xylar Asay-Davis
fcMOC = moc(gf)
if geojson_filename is not None:
fcMOC.to_geojson(geojson_filename)
dsMesh = xarray.open_dataset(mesh_filename)
dsMasks = mpas_tools.mesh.conversion.mask(
dsMesh=dsMesh, fcMask=fcMOC, logger=logger, dir=dir
)
if mask_filename is not None:
write_netcdf(dsMasks, mask_filename, char_dim_name='StrLen')
dsMasksAndTransects = add_moc_southern_boundary_transects(
dsMasks, dsMesh, logger=logger
)
write_netcdf(
dsMasksAndTransects, mask_and_transect_filename, char_dim_name='StrLen'
)
[docs]
def add_moc_southern_boundary_transects(dsMask, dsMesh, logger=None):
"""
Parameters
----------
dsMask : ``xarray.Dataset``
Region masks defining MOC basins
dsMesh : ``xarray.Dataset``, optional
An MPAS mesh on which the masks should be created
logger : ``logging.Logger``, optional
A logger for the output if not stdout
Returns
-------
dsMask : ``xarray.Dataset``
Region masks defining MOC basins and the corresponding
southern-boundary transects
"""
useStdout = logger is None
if useStdout:
logger = logging.getLogger()
logger.addHandler(logging.StreamHandler(sys.stdout))
logger.setLevel(logging.INFO)
(
southernBoundaryEdges,
southernBoundaryEdgeSigns,
southernBoundaryVertices,
) = _extract_southern_boundary(
dsMesh, dsMask, latBuffer=3.0 * numpy.pi / 180.0, logger=logger
)
_add_transects_to_moc(
dsMesh,
dsMask,
southernBoundaryEdges,
southernBoundaryEdgeSigns,
southernBoundaryVertices,
)
if useStdout:
logger.handlers = []
return dsMask
def moc_southern_boundary_extractor():
"""
Entry point for moc_southern_boundary_extractor command-line tool
"""
description = """
This script takes a mesh file (-m flag) and a file with MOC regions masks
(-f flag) produce by the MPAS mask creator. The script produces a copy of
the contents of the MOC mask file, adding transects that mark the southern
boundary of each region in a file indicated with the -o flag. The transect
is applied only to vertices and edges, not cells, because the need for southern
boundary transect data on cells is not foreseen.
"""
parser = argparse.ArgumentParser(
description=description, formatter_class=argparse.RawTextHelpFormatter
)
parser.add_argument(
'-f',
'--in_file',
dest='in_file',
help='Input file with MOC masks',
metavar='IN_FILE',
required=True,
)
parser.add_argument(
'-m',
'--mesh_file',
dest='mesh_file',
help='Input mesh file',
metavar='MESH_FILE',
required=True,
)
parser.add_argument(
'-o',
'--out_file',
dest='out_file',
help='Output file for MOC masks and southern-boundary transects',
metavar='OUT_FILE',
required=True,
)
args = parser.parse_args()
dsMasks = xarray.open_dataset(args.in_file)
dsMesh = xarray.open_dataset(args.mesh_file)
dsMasksAndTransects = add_moc_southern_boundary_transects(dsMasks, dsMesh)
write_netcdf(dsMasksAndTransects, args.out_file)
def _extract_southern_boundary(mesh, mocMask, latBuffer, logger):
"""
Extracts the southern boundary of each region mask in mocMask. Mesh info
is taken from mesh. latBuffer is a number of radians above the southern-
most point that should be considered to definitely be in the southern
boundary.
"""
nCells = mesh.sizes['nCells']
nEdges = mesh.sizes['nEdges']
nRegions = mocMask.sizes['nRegions']
assert mocMask.sizes['nCells'] == nCells
# convert to python zero-based indices
cellsOnEdge = mesh.variables['cellsOnEdge'].values - 1
verticesOnEdge = mesh.variables['verticesOnEdge'].values - 1
edgesOnVertex = mesh.variables['edgesOnVertex'].values - 1
latEdge = mesh.variables['latEdge'].values
cellsOnEdgeInRange = numpy.logical_and(
cellsOnEdge >= 0, cellsOnEdge < nCells
)
southernBoundaryEdges = []
southernBoundaryEdgeSigns = []
southernBoundaryVertices = []
for iRegion in range(nRegions):
name = mocMask.regionNames[iRegion].values.astype('U')
logger.info(name)
cellMask = mocMask.variables['regionCellMasks'][:, iRegion].values
# land cells are outside not in the MOC region
cellsOnEdgeMask = numpy.zeros(cellsOnEdge.shape, bool)
# set mask values for cells that are in range (not land)
cellsOnEdgeMask[cellsOnEdgeInRange] = (
cellMask[cellsOnEdge[cellsOnEdgeInRange]] == 1
)
logger.info(' computing edge sign...')
edgeSign = numpy.zeros(nEdges)
# positive sign if the first cell on edge is in the region
mask = numpy.logical_and(
cellsOnEdgeMask[:, 0], numpy.logical_not(cellsOnEdgeMask[:, 1])
)
edgeSign[mask] = -1.0
# negative sign if the second cell on edge is in the region
mask = numpy.logical_and(
cellsOnEdgeMask[:, 1], numpy.logical_not(cellsOnEdgeMask[:, 0])
)
edgeSign[mask] = 1.0
isMOCBoundaryEdge = edgeSign != 0.0
edgesMOCBoundary = numpy.arange(nEdges)[isMOCBoundaryEdge]
logger.info(' done.')
startEdge = numpy.argmin(latEdge[isMOCBoundaryEdge])
startEdge = edgesMOCBoundary[startEdge]
minLat = latEdge[startEdge]
logger.info(' getting edge sequence...')
# follow the boundary from this point to get a loop of edges
# Note: it is possible but unlikely that the southern-most point is
# not within bulk region of the MOC mask if the region is not a single
# shape
edgeSequence, edgeSequenceSigns, vertexSequence = (
_get_edge_sequence_on_boundary(
startEdge, edgeSign, edgesOnVertex, verticesOnEdge
)
)
logger.info(f' done: {len(edgeSequence)} edges in transect.')
aboveSouthernBoundary = latEdge[edgeSequence] > minLat + latBuffer
# find and eliminate the longest contiguous sequence (if any) from the
# edge sequence that is above the possible region of the soutehrn
# boundary
belowToAbove = numpy.logical_and(
numpy.logical_not(aboveSouthernBoundary[0:-1]),
aboveSouthernBoundary[1:],
)
aboveToBelow = numpy.logical_and(
aboveSouthernBoundary[0:-1],
numpy.logical_not(aboveSouthernBoundary[1:]),
)
startIndices = numpy.arange(1, len(edgeSequence))[belowToAbove]
endIndices = numpy.arange(1, len(edgeSequence))[aboveToBelow]
assert len(startIndices) == len(endIndices)
if len(startIndices) == 0:
# the whole sequence is the southern boundary
southernBoundaryEdges.append(edgeSequence)
southernBoundaryEdgeSigns.append(edgeSequenceSigns)
southernBoundaryVertices.append(vertexSequence)
continue
# there are some parts of the sequence above the boundary. Let's
# eliminate the longest one.
aboveLength = endIndices - startIndices
longest = numpy.argmax(aboveLength)
# we want all the indices in the sequence *not* part of the longest
indices = numpy.arange(
endIndices[longest], startIndices[longest] + len(edgeSequence)
)
indices = numpy.mod(indices, len(edgeSequence))
southernBoundaryEdges.append(edgeSequence[indices])
southernBoundaryEdgeSigns.append(edgeSequenceSigns[indices])
# we want one extra vertex in the vertex sequence
indices = numpy.arange(
endIndices[longest], startIndices[longest] + len(edgeSequence) + 1
)
indices = numpy.mod(indices, len(edgeSequence))
southernBoundaryVertices.append(vertexSequence[indices])
return (
southernBoundaryEdges,
southernBoundaryEdgeSigns,
southernBoundaryVertices,
)
def _add_transects_to_moc(
mesh,
mocMask,
southernBoundaryEdges,
southernBoiundaryEdgeSigns,
southernBoundaryVertices,
):
"""
Creates transect fields in mocMask from the edges, edge signs and
vertices defining the southern boundaries. Mesh info (nEdges and
nVertices) is taken from the mesh file.
"""
nTransects = len(southernBoundaryEdges)
nEdges = mesh.sizes['nEdges']
nVertices = mesh.sizes['nVertices']
maxEdgesInTransect = numpy.amax(
[
len(southernBoundaryEdges[iTransect])
for iTransect in range(nTransects)
]
)
maxVerticesInTransect = numpy.amax(
[
len(southernBoundaryVertices[iTransect])
for iTransect in range(nTransects)
]
)
transectEdgeMasks = numpy.zeros((nEdges, nTransects), numpy.int32)
transectEdgeMaskSigns = numpy.zeros((nEdges, nTransects), numpy.int32)
transectEdgeGlobalIDs = numpy.zeros(
(nTransects, maxEdgesInTransect), numpy.int32
)
transectVertexMasks = numpy.zeros((nVertices, nTransects), numpy.int32)
transectVertexGlobalIDs = numpy.zeros(
(nTransects, maxVerticesInTransect), numpy.int32
)
for iTransect in range(nTransects):
transectEdgeMasks[southernBoundaryEdges[iTransect], iTransect] = 1
transectEdgeMaskSigns[southernBoundaryEdges[iTransect], iTransect] = (
southernBoiundaryEdgeSigns[iTransect]
)
transectCount = len(southernBoundaryEdges[iTransect])
transectEdgeGlobalIDs[iTransect, 0:transectCount] = (
southernBoundaryEdges[iTransect] + 1
)
transectVertexMasks[southernBoundaryVertices[iTransect], iTransect] = 1
transectCount = len(southernBoundaryVertices[iTransect])
transectVertexGlobalIDs[iTransect, 0:transectCount] = (
southernBoundaryVertices[iTransect] + 1
)
mocMask['transectEdgeMasks'] = (
('nEdges', 'nTransects'),
transectEdgeMasks,
)
mocMask['transectEdgeMaskSigns'] = (
('nEdges', 'nTransects'),
transectEdgeMaskSigns,
)
mocMask['transectEdgeGlobalIDs'] = (
('nTransects', 'maxEdgesInTransect'),
transectEdgeGlobalIDs,
)
mocMask['transectVertexMasks'] = (
('nVertices', 'nTransects'),
transectVertexMasks,
)
mocMask['transectVertexGlobalIDs'] = (
('nTransects', 'maxVerticesInTransect'),
transectVertexGlobalIDs,
)
if 'nRegionsInGroup' not in mocMask:
nRegions = mocMask.sizes['nRegions']
nRegionGroups = 2
nRegionsInGroup = nRegions * numpy.ones(
nRegionGroups, dtype=numpy.int32
)
regionsInGroup = numpy.zeros(
(nRegionGroups, nRegions), dtype=numpy.int32
)
regionGroupNames = ['MOCBasinRegionsGroup', 'all']
regionNames = mocMask.regionNames.values
for index in range(nRegionGroups):
regionsInGroup[index, :] = numpy.arange(1, nRegions + 1)
mocMask['nRegionsInGroup'] = (('nRegionGroups',), nRegionsInGroup)
mocMask['regionsInGroup'] = (
('nRegionGroups', 'maxRegionsInGroup'),
regionsInGroup,
)
mocMask['regionGroupNames'] = (
('nRegionGroups',),
numpy.zeros((nRegionGroups,), dtype=f'|S{default_nchar}'),
)
for index in range(nRegionGroups):
mocMask['regionGroupNames'][index] = regionGroupNames[index]
# we need to make sure the region names use the same string length
mocMask['regionNames'] = (
('nRegions',),
numpy.zeros((nRegions,), dtype=f'|S{default_nchar}'),
)
for index in range(nRegions):
mocMask['regionNames'][index] = regionNames[index]
mocMask['transectNames'] = mocMask.regionNames.rename(
{'nRegions': 'nTransects'}
)
mocMask['nTransectsInGroup'] = mocMask.nRegionsInGroup.rename(
{'nRegionGroups': 'nTransectGroups'}
)
mocMask['transectsInGroup'] = mocMask.regionsInGroup.rename(
{
'nRegionGroups': 'nTransectGroups',
'maxRegionsInGroup': 'maxTransectsInGroup',
}
)
mocMask['transectGroupNames'] = mocMask.regionGroupNames.rename(
{'nRegionGroups': 'nTransectGroups'}
)
def _get_edge_sequence_on_boundary(
startEdge, edgeSign, edgesOnVertex, verticesOnEdge
):
"""
Follows the boundary from a starting edge to produce a sequence of
edges that form a closed loop.
startEdge is an edge on the boundary that will be both the start and
end of the loop.
isBoundaryEdge is a mask that indicates which edges are on the
boundary
returns lists of edges, edge signs and vertices
"""
iEdge = startEdge
edgeSequence = []
vertexSequence = []
while True:
assert edgeSign[iEdge] == 1.0 or edgeSign[iEdge] == -1.0
if edgeSign[iEdge] == 1.0:
v = 0
else:
v = 1
iVertex = verticesOnEdge[iEdge, v]
eov = edgesOnVertex[iVertex, :]
# find the edge that is not iEdge but is on the boundary
nextEdge = -1
for edge in eov:
if edge != -1 and edge != iEdge and edgeSign[edge] != 0:
nextEdge = edge
break
assert nextEdge != -1
edgeSequence.append(iEdge)
vertexSequence.append(iVertex)
iEdge = nextEdge
if iEdge == startEdge:
break
edgeSequence = numpy.array(edgeSequence)
edgeSequenceSigns = edgeSign[edgeSequence]
vertexSequence = numpy.array(vertexSequence)
return edgeSequence, edgeSequenceSigns, vertexSequence