From ec1d07cc8310903a1cea25723b8df2bd974b2d0f Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Wed, 5 Feb 2025 13:22:17 -0600 Subject: [PATCH 1/3] Make sure we use int32 in mask creation --- conda_package/mpas_tools/mesh/mask.py | 28 +++++++++++++-------------- 1 file changed, 14 insertions(+), 14 deletions(-) diff --git a/conda_package/mpas_tools/mesh/mask.py b/conda_package/mpas_tools/mesh/mask.py index c14e20e37..3a91bf6d4 100644 --- a/conda_package/mpas_tools/mesh/mask.py +++ b/conda_package/mpas_tools/mesh/mask.py @@ -100,11 +100,11 @@ 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 @@ -276,12 +276,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] @@ -290,7 +290,7 @@ 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]) @@ -447,7 +447,7 @@ 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.') @@ -559,12 +559,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',), @@ -719,12 +719,12 @@ 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',), @@ -1180,7 +1180,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): @@ -1248,7 +1248,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 = [] @@ -1272,7 +1272,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)): @@ -1300,7 +1300,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 From b937f16f88d415797a90afc28076817a85ac3b7a Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Wed, 5 Feb 2025 13:29:54 -0600 Subject: [PATCH 2/3] Lint the code --- conda_package/mpas_tools/mesh/mask.py | 54 +++++++++++++++------------ 1 file changed, 31 insertions(+), 23 deletions(-) diff --git a/conda_package/mpas_tools/mesh/mask.py b/conda_package/mpas_tools/mesh/mask.py index 3a91bf6d4..ebb016495 100644 --- a/conda_package/mpas_tools/mesh/mask.py +++ b/conda_package/mpas_tools/mesh/mask.py @@ -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 ---------- @@ -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=numpy.int32)) + ((dim, 'nRegions'), numpy.zeros((nPoints, nRegions), + dtype=numpy.int32)) for index in range(nRegions): mask = masks[index] - dsMasks[masksVarName][:, index] = numpy.array(mask, dtype=numpy.int32) + dsMasks[masksVarName][:, index] = numpy.array(mask, + dtype=numpy.int32) if 'regionNames' not in dsMasks: # create a new data array for mask names @@ -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=numpy.int32) + dsMasks[masksVarName][:, index] = numpy.array(mask, + dtype=numpy.int32) if addEdgeSign and maskType == 'edge': print(transectNames[index]) @@ -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=numpy.int32)) + dsMasks[masksVarName] = (('nCells',), numpy.array(seedMask, + dtype=numpy.int32)) if logger is not None: logger.info(' Done.') @@ -719,7 +726,8 @@ 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=numpy.int32)) + ((ydim, xdim, 'nRegions'), numpy.zeros((ny, nx, nRegions), + dtype=numpy.int32)) for index in range(nRegions): mask = masks[index] @@ -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 @@ -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] From 8e3eade699ebbcc43c4de4befd48a9677a67a836 Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Wed, 5 Feb 2025 13:33:22 -0600 Subject: [PATCH 3/3] Make sure to use int32 in MOC southern transect --- conda_package/mpas_tools/ocean/moc.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/conda_package/mpas_tools/ocean/moc.py b/conda_package/mpas_tools/ocean/moc.py index 16d2d45e1..afa8544e5 100755 --- a/conda_package/mpas_tools/ocean/moc.py +++ b/conda_package/mpas_tools/ocean/moc.py @@ -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 @@ -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