Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Make sure we use int32 in mask creation #603

Open
wants to merge 3 commits into
base: master
Choose a base branch
from
Open

Make sure we use int32 in mask creation #603

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
ec1d07c
Make sure we use int32 in mask creation
xylar Feb 5, 2025
b937f16
Lint the code
xylar Feb 5, 2025
8e3eade
Make sure to use int32 in MOC southern transect
xylar Feb 5, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
72 changes: 40 additions & 32 deletions conda_package/mpas_tools/mesh/mask.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,31 +1,34 @@
import xarray as xr
import argparse
from functools import partial

import numpy
from scipy.spatial import KDTree
import progressbar
import shapely.geometry
import shapely.ops
from shapely.geometry import box, Polygon, MultiPolygon, GeometryCollection
from shapely.strtree import STRtree
import progressbar
from functools import partial
import argparse
import xarray as xr
from igraph import Graph
from scipy.spatial import KDTree
from shapely.geometry import GeometryCollection, MultiPolygon, Polygon, box
from shapely.strtree import STRtree

from geometric_features import read_feature_collection

from mpas_tools.transects import subdivide_great_circle, \
lon_lat_to_cartesian, cartesian_to_lon_lat
from mpas_tools.parallel import create_pool
from mpas_tools.cime.constants import constants
from mpas_tools.io import write_netcdf
from mpas_tools.logging import LoggingContext
from mpas_tools.cime.constants import constants
from mpas_tools.parallel import create_pool
from mpas_tools.transects import (
cartesian_to_lon_lat,
lon_lat_to_cartesian,
subdivide_great_circle,
)


def compute_mpas_region_masks(dsMesh, fcMask, maskTypes=('cell', 'vertex'),
logger=None, pool=None, chunkSize=1000,
showProgress=False, subdivisionThreshold=30.):
"""
Use shapely and processes to create a set of masks from a feature collection
made up of regions (polygons)
Use shapely and processes to create a set of masks from a feature
collection made up of regions (polygons)

Parameters
----------
Expand Down Expand Up @@ -100,11 +103,13 @@ def compute_mpas_region_masks(dsMesh, fcMask, maskTypes=('cell', 'vertex'),
# create a new data array for masks
masksVarName = 'region{}Masks'.format(suffix)
dsMasks[masksVarName] = \
((dim, 'nRegions'), numpy.zeros((nPoints, nRegions), dtype=int))
((dim, 'nRegions'), numpy.zeros((nPoints, nRegions),
dtype=numpy.int32))

for index in range(nRegions):
mask = masks[index]
dsMasks[masksVarName][:, index] = numpy.array(mask, dtype=int)
dsMasks[masksVarName][:, index] = numpy.array(mask,
dtype=numpy.int32)

if 'regionNames' not in dsMasks:
# create a new data array for mask names
Expand Down Expand Up @@ -276,12 +281,12 @@ def compute_mpas_transect_masks(dsMesh, fcMask, earthRadius,
masksVarName = 'transect{}Masks'.format(suffix)
dsMasks[masksVarName] = \
((dim, 'nTransects'),
numpy.zeros((nPolygons, nTransects), dtype=int))
numpy.zeros((nPolygons, nTransects), dtype=numpy.int32))

if addEdgeSign and maskType == 'edge':
dsMasks['transectEdgeMaskSigns'] = \
((dim, 'nTransects'),
numpy.zeros((nPolygons, nTransects), dtype=int))
numpy.zeros((nPolygons, nTransects), dtype=numpy.int32))

for index in range(nTransects):
maskAndDuplicates = masks[index]
Expand All @@ -290,7 +295,8 @@ def compute_mpas_transect_masks(dsMesh, fcMask, earthRadius,
mask[duplicatePolygons] = \
numpy.logical_or(mask[duplicatePolygons],
maskAndDuplicates[nPolygons:])
dsMasks[masksVarName][:, index] = numpy.array(mask, dtype=int)
dsMasks[masksVarName][:, index] = numpy.array(mask,
dtype=numpy.int32)

if addEdgeSign and maskType == 'edge':
print(transectNames[index])
Expand Down Expand Up @@ -447,7 +453,8 @@ def compute_mpas_flood_fill_mask(dsMesh, fcSeed, daGrow=None, logger=None,
logger.info(' Adding masks to dataset...')
# create a new data array for the mask
masksVarName = 'cellSeedMask'
dsMasks[masksVarName] = (('nCells',), numpy.array(seedMask, dtype=int))
dsMasks[masksVarName] = (('nCells',), numpy.array(seedMask,
dtype=numpy.int32))

if logger is not None:
logger.info(' Done.')
Expand Down Expand Up @@ -559,12 +566,12 @@ def compute_lon_lat_region_masks(lon, lat, fcMask, logger=None, pool=None,
masksVarName = 'regionMasks'
dsMasks[masksVarName] = \
(('lat', 'lon', 'nRegions'), numpy.zeros((nlat, nlon, nRegions),
dtype=int))
dtype=numpy.int32))

for index in range(nRegions):
mask = masks[index]
dsMasks[masksVarName][:, :, index] = \
numpy.array(mask.reshape(shape), dtype=int)
numpy.array(mask.reshape(shape), dtype=numpy.int32)

# create a new data array for mask names
dsMasks['regionNames'] = (('nRegions',),
Expand Down Expand Up @@ -719,12 +726,13 @@ def compute_projection_grid_region_masks(
# create a new data array for masks
masksVarName = 'regionMasks'
dsMasks[masksVarName] = \
((ydim, xdim, 'nRegions'), numpy.zeros((ny, nx, nRegions), dtype=int))
((ydim, xdim, 'nRegions'), numpy.zeros((ny, nx, nRegions),
dtype=numpy.int32))

for index in range(nRegions):
mask = masks[index]
dsMasks[masksVarName][:, :, index] = \
numpy.array(mask.reshape((ny, nx)), dtype=int)
numpy.array(mask.reshape((ny, nx)), dtype=numpy.int32)

# create a new data array for mask names
dsMasks['regionNames'] = (('nRegions',),
Expand Down Expand Up @@ -932,9 +940,9 @@ def _katana(geometry, threshold, count=0, maxcount=250):

1. Redistributions of source code must retain the above copyright notice,
this list of conditions and the following disclaimer.
2. Redistributions in binary form must reproduce the above copyright notice,
this list of conditions and the following disclaimer in the documentation
and/or other materials provided with the distribution.
2. Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in the
documentation and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS “AS IS”
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
Expand Down Expand Up @@ -1074,7 +1082,7 @@ def _get_polygons(dsMesh, maskType):
for iVertex in range(1, maxVertices):
mask = firstValid < 0
firstValid[mask] = vertexIndices[mask, iVertex]
assert(numpy.all(firstValid >= 0))
assert numpy.all(firstValid >= 0)
for iVertex in range(maxVertices):
mask = vertexIndices[:, iVertex] < 0
vertexIndices[mask, iVertex] = firstValid[mask]
Expand Down Expand Up @@ -1180,7 +1188,7 @@ def _compute_seed_mask(fcSeed, lon, lat, workers):

tree = KDTree(points)

mask = numpy.zeros(len(lon), dtype=int)
mask = numpy.zeros(len(lon), dtype=numpy.int32)

points = numpy.zeros((len(fcSeed.features), 2))
for index, feature in enumerate(fcSeed.features):
Expand Down Expand Up @@ -1248,7 +1256,7 @@ def _compute_edge_sign(dsMesh, edgeMask, shape):

unique_vertices = numpy.unique(voe.ravel())

local_voe = numpy.zeros(voe.shape, dtype=int)
local_voe = numpy.zeros(voe.shape, dtype=numpy.int32)
distance = []
unique_lon = []
unique_lat = []
Expand All @@ -1272,7 +1280,7 @@ def _compute_edge_sign(dsMesh, edgeMask, shape):
graph.es['edges'] = edge_indices
graph.es['vertices'] = [(v0, v1) for v0, v1 in zip(voe[:, 0], voe[:, 1])]

edgeSign = numpy.zeros(edgeMask.shape, dtype=int)
edgeSign = numpy.zeros(edgeMask.shape, dtype=numpy.int32)

clusters = graph.connected_components()
for cluster_index in range(len(clusters)):
Expand Down Expand Up @@ -1300,7 +1308,7 @@ def _compute_edge_sign(dsMesh, edgeMask, shape):
sign = [-1]
else:
verts = numpy.array(edges['vertices'])
sign = numpy.zeros(len(indices), dtype=int)
sign = numpy.zeros(len(indices), dtype=numpy.int32)
for index in range(len(indices)-1):
if verts[index, 1] in verts[index+1, :]:
sign[index] = 1
Expand Down
7 changes: 4 additions & 3 deletions conda_package/mpas_tools/ocean/moc.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,8 +2,8 @@
import logging
import sys

import xarray
import numpy
import xarray
from geometric_features.aggregation.ocean import moc

import mpas_tools.mesh.conversion
Expand Down Expand Up @@ -328,8 +328,9 @@ def _add_transects_to_moc(mesh, mocMask, southernBoundaryEdges,
if 'nRegionsInGroup' not in mocMask:
nRegions = mocMask.sizes['nRegions']
nRegionGroups = 2
nRegionsInGroup = nRegions*numpy.ones(nRegionGroups, dtype=int)
regionsInGroup = numpy.zeros((nRegionGroups, nRegions), dtype=int)
nRegionsInGroup = nRegions*numpy.ones(nRegionGroups, dtype=numpy.int32)
regionsInGroup = numpy.zeros((nRegionGroups, nRegions),
dtype=numpy.int32)
regionGroupNames = ['MOCBasinRegionsGroup', 'all']
regionNames = mocMask.regionNames.values
nChar = 64
Expand Down
Loading