From 82927be18be9cf3c64f56f86e4c2b3f43ac18479 Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Sat, 16 Sep 2023 16:49:11 +0200 Subject: [PATCH 1/6] Add MpasCellMeshDescriptor Deprecate MpasMeshDescriptor --- pyremap/descriptor/__init__.py | 1 + .../descriptor/mpas_cell_mesh_descriptor.py | 203 ++++++++++++++++++ pyremap/descriptor/mpas_mesh_descriptor.py | 166 +------------- 3 files changed, 211 insertions(+), 159 deletions(-) create mode 100644 pyremap/descriptor/mpas_cell_mesh_descriptor.py diff --git a/pyremap/descriptor/__init__.py b/pyremap/descriptor/__init__.py index b30c3fd..df40e1b 100644 --- a/pyremap/descriptor/__init__.py +++ b/pyremap/descriptor/__init__.py @@ -17,6 +17,7 @@ get_lat_lon_descriptor, ) from pyremap.descriptor.mesh_descriptor import MeshDescriptor +from pyremap.descriptor.mpas_cell_mesh_descriptor import MpasCellMeshDescriptor from pyremap.descriptor.mpas_edge_mesh_descriptor import MpasEdgeMeshDescriptor from pyremap.descriptor.mpas_mesh_descriptor import MpasMeshDescriptor from pyremap.descriptor.point_collection_descriptor import ( diff --git a/pyremap/descriptor/mpas_cell_mesh_descriptor.py b/pyremap/descriptor/mpas_cell_mesh_descriptor.py new file mode 100644 index 0000000..71fac7b --- /dev/null +++ b/pyremap/descriptor/mpas_cell_mesh_descriptor.py @@ -0,0 +1,203 @@ +# This software is open source software available under the BSD-3 license. +# +# Copyright (c) 2022 Triad National Security, LLC. All rights reserved. +# Copyright (c) 2022 Lawrence Livermore National Security, LLC. All rights +# reserved. +# Copyright (c) 2022 UT-Battelle, LLC. All rights reserved. +# +# Additional copyright and license information can be found in the LICENSE file +# distributed with this code, or at +# https://raw.githubusercontent.com/MPAS-Dev/pyremap/main/LICENSE + +import sys +import warnings + +import netCDF4 +import numpy +import xarray + +from pyremap.descriptor.mesh_descriptor import MeshDescriptor +from pyremap.descriptor.utility import create_scrip + + +class MpasCellMeshDescriptor(MeshDescriptor): + """ + A class for describing an MPAS cell mesh + + Attributes + ---------- + vertices : bool + Whether the mapping is to or from vertices instead of corners + (for non-conservative remapping) + + fileName : str + The path of the file containing the MPAS mesh + """ + def __init__(self, fileName, meshName=None, vertices=False): + """ + Constructor stores the file name + + Parameters + ---------- + fileName : str + The path of the file containing the MPAS mesh + + meshName : str, optional + The name of the MPAS mesh (e.g. ``'oEC60to30'`` or + ``'oRRS18to6'``). If not provided, the data set in ``fileName`` + must have a global attribute ``meshName`` that will be used + instead. + + vertices : bool, optional + Whether the mapping is to or from vertices instead of corners + (for non-conservative remapping) + """ + super().__init__() + + if vertices: + warnings.warn('Creating a MpasCellMeshDescriptor with ' + 'vertices=True is deprecated and will be removed in ' + 'the next release. Use MpasVertexMeshDescriptor ' + 'instead.', DeprecationWarning) + + self.vertices = vertices + + with xarray.open_dataset(fileName) as ds: + + if meshName is None: + if 'meshName' not in ds.attrs: + raise ValueError('No meshName provided or found in file.') + self.meshName = ds.attrs['meshName'] + else: + self.meshName = meshName + + self.fileName = fileName + self.regional = True + + if vertices: + # build coords + self.coords = {'latVertex': {'dims': 'nVertices', + 'data': ds.latVertex.values, + 'attrs': {'units': 'radians'}}, + 'lonVertex': {'dims': 'nVertices', + 'data': ds.lonVertex.values, + 'attrs': {'units': 'radians'}}} + self.dims = ['nVertices'] + else: + # build coords + self.coords = {'latCell': {'dims': 'nCells', + 'data': ds.latCell.values, + 'attrs': {'units': 'radians'}}, + 'lonCell': {'dims': 'nCells', + 'data': ds.lonCell.values, + 'attrs': {'units': 'radians'}}} + self.dims = ['nCells'] + self.dimSize = [ds.dims[dim] for dim in self.dims] + + def to_scrip(self, scripFileName): + """ + Given an MPAS mesh file, create a SCRIP file based on the mesh. + + Parameters + ---------- + scripFileName : str + The path to which the SCRIP file should be written + """ + if self.vertices: + raise ValueError('A SCRIP file won\'t work for remapping vertices') + + inFile = netCDF4.Dataset(self.fileName, 'r') + outFile = netCDF4.Dataset(scripFileName, 'w', format=self.format) + + # Get info from input file + latCell = inFile.variables['latCell'][:] + lonCell = inFile.variables['lonCell'][:] + latVertex = inFile.variables['latVertex'][:] + lonVertex = inFile.variables['lonVertex'][:] + verticesOnCell = inFile.variables['verticesOnCell'][:] + nEdgesOnCell = inFile.variables['nEdgesOnCell'][:] + nCells = len(inFile.dimensions['nCells']) + maxVertices = len(inFile.dimensions['maxEdges']) + areaCell = inFile.variables['areaCell'][:] + sphereRadius = float(inFile.sphere_radius) + + create_scrip(outFile, grid_size=nCells, grid_corners=maxVertices, + grid_rank=1, units='radians', meshName=self.meshName) + + grid_area = outFile.createVariable('grid_area', 'f8', ('grid_size',)) + grid_area.units = 'radian^2' + # SCRIP uses square radians + grid_area[:] = areaCell[:] / (sphereRadius**2) + + outFile.variables['grid_center_lat'][:] = latCell[:] + outFile.variables['grid_center_lon'][:] = lonCell[:] + outFile.variables['grid_dims'][:] = nCells + outFile.variables['grid_imask'][:] = 1 + + # grid corners: + grid_corner_lon = numpy.zeros((nCells, maxVertices)) + grid_corner_lat = numpy.zeros((nCells, maxVertices)) + for iVertex in range(maxVertices): + cellIndices = numpy.arange(nCells) + # repeat the last vertex wherever iVertex > nEdgesOnCell + localVertexIndices = numpy.minimum(nEdgesOnCell - 1, iVertex) + vertexIndices = verticesOnCell[cellIndices, localVertexIndices] - 1 + grid_corner_lat[cellIndices, iVertex] = latVertex[vertexIndices] + grid_corner_lon[cellIndices, iVertex] = lonVertex[vertexIndices] + + outFile.variables['grid_corner_lat'][:] = grid_corner_lat[:] + outFile.variables['grid_corner_lon'][:] = grid_corner_lon[:] + + # Update history attribute of netCDF file + if hasattr(inFile, 'history'): + newhist = '\n'.join([getattr(inFile, 'history'), + ' '.join(sys.argv[:])]) + else: + newhist = sys.argv[:] + setattr(outFile, 'history', newhist) + + inFile.close() + outFile.close() + + def to_esmf(self, esmfFileName): + """ + Create an ESMF mesh file for the mesh + + Parameters + ---------- + esmfFileName : str + The path to which the ESMF mesh file should be written + """ + with xarray.open_dataset(self.fileName) as ds: + + nodeCount = ds.sizes['nVertices'] + elementCount = ds.sizes['nCells'] + coordDim = 2 + + nodeCoords = numpy.zeros((nodeCount, coordDim)) + nodeCoords[:, 0] = ds.lonVertex.values + nodeCoords[:, 1] = ds.latVertex.values + centerCoords = numpy.zeros((elementCount, coordDim)) + centerCoords[:, 0] = ds.lonCell.values + centerCoords[:, 1] = ds.latCell.values + + elementConn = ds.verticesOnCell.values + elementConn[elementConn == nodeCount + 1] = -1 + + ds_out = xarray.Dataset() + ds_out['nodeCoords'] = (('nodeCount', 'coordDim'), nodeCoords) + ds_out.nodeCoords.attrs['units'] = 'radians' + ds_out['centerCoords'] = \ + (('elementCount', 'coordDim'), centerCoords) + ds_out.centerCoords.attrs['units'] = 'radians' + ds_out['elementConn'] = \ + (('elementCount', 'maxNodePElement'), elementConn) + ds_out.elementConn.attrs['start_index'] = 1 + ds_out.elementConn.attrs['_FillValue'] = -1 + ds_out['numElementConn'] = \ + (('elementCount',), ds.nEdgesOnCell.values) + ds_out.attrs['gridType'] = 'unstructured mesh' + ds_out.attrs['version'] = '0.9' + + ds_out.to_netcdf(esmfFileName, format=self.format, + engine=self.engine) diff --git a/pyremap/descriptor/mpas_mesh_descriptor.py b/pyremap/descriptor/mpas_mesh_descriptor.py index 2204a18..f8b0671 100644 --- a/pyremap/descriptor/mpas_mesh_descriptor.py +++ b/pyremap/descriptor/mpas_mesh_descriptor.py @@ -9,29 +9,15 @@ # distributed with this code, or at # https://raw.githubusercontent.com/MPAS-Dev/pyremap/main/LICENSE -import sys +import warnings -import netCDF4 -import numpy -import xarray +from pyremap.descriptor.mpas_cell_mesh_descriptor import MpasCellMeshDescriptor -from pyremap.descriptor.mesh_descriptor import MeshDescriptor -from pyremap.descriptor.utility import create_scrip - -class MpasMeshDescriptor(MeshDescriptor): +class MpasMeshDescriptor(MpasCellMeshDescriptor): """ A class for describing an MPAS cell mesh (which also supports non-conservative remapping from vertices) - - Attributes - ---------- - vertices : bool - Whether the mapping is to or from vertices instead of corners - (for non-conservative remapping) - - fileName : str - The path of the file containing the MPAS mesh """ def __init__(self, fileName, meshName=None, vertices=False): """ @@ -52,146 +38,8 @@ def __init__(self, fileName, meshName=None, vertices=False): Whether the mapping is to or from vertices instead of corners (for non-conservative remapping) """ - super().__init__() - - self.vertices = vertices - - with xarray.open_dataset(fileName) as ds: - - if meshName is None: - if 'meshName' not in ds.attrs: - raise ValueError('No meshName provided or found in file.') - self.meshName = ds.attrs['meshName'] - else: - self.meshName = meshName - - self.fileName = fileName - self.regional = True - - if vertices: - # build coords - self.coords = {'latVertex': {'dims': 'nVertices', - 'data': ds.latVertex.values, - 'attrs': {'units': 'radians'}}, - 'lonVertex': {'dims': 'nVertices', - 'data': ds.lonVertex.values, - 'attrs': {'units': 'radians'}}} - self.dims = ['nVertices'] - else: - # build coords - self.coords = {'latCell': {'dims': 'nCells', - 'data': ds.latCell.values, - 'attrs': {'units': 'radians'}}, - 'lonCell': {'dims': 'nCells', - 'data': ds.lonCell.values, - 'attrs': {'units': 'radians'}}} - self.dims = ['nCells'] - self.dimSize = [ds.dims[dim] for dim in self.dims] - - def to_scrip(self, scripFileName): - """ - Given an MPAS mesh file, create a SCRIP file based on the mesh. - - Parameters - ---------- - scripFileName : str - The path to which the SCRIP file should be written - """ - if self.vertices: - raise ValueError('A SCRIP file won\'t work for remapping vertices') - - inFile = netCDF4.Dataset(self.fileName, 'r') - outFile = netCDF4.Dataset(scripFileName, 'w', format=self.format) - - # Get info from input file - latCell = inFile.variables['latCell'][:] - lonCell = inFile.variables['lonCell'][:] - latVertex = inFile.variables['latVertex'][:] - lonVertex = inFile.variables['lonVertex'][:] - verticesOnCell = inFile.variables['verticesOnCell'][:] - nEdgesOnCell = inFile.variables['nEdgesOnCell'][:] - nCells = len(inFile.dimensions['nCells']) - maxVertices = len(inFile.dimensions['maxEdges']) - areaCell = inFile.variables['areaCell'][:] - sphereRadius = float(inFile.sphere_radius) - - create_scrip(outFile, grid_size=nCells, grid_corners=maxVertices, - grid_rank=1, units='radians', meshName=self.meshName) - - grid_area = outFile.createVariable('grid_area', 'f8', ('grid_size',)) - grid_area.units = 'radian^2' - # SCRIP uses square radians - grid_area[:] = areaCell[:] / (sphereRadius**2) - - outFile.variables['grid_center_lat'][:] = latCell[:] - outFile.variables['grid_center_lon'][:] = lonCell[:] - outFile.variables['grid_dims'][:] = nCells - outFile.variables['grid_imask'][:] = 1 - - # grid corners: - grid_corner_lon = numpy.zeros((nCells, maxVertices)) - grid_corner_lat = numpy.zeros((nCells, maxVertices)) - for iVertex in range(maxVertices): - cellIndices = numpy.arange(nCells) - # repeat the last vertex wherever iVertex > nEdgesOnCell - localVertexIndices = numpy.minimum(nEdgesOnCell - 1, iVertex) - vertexIndices = verticesOnCell[cellIndices, localVertexIndices] - 1 - grid_corner_lat[cellIndices, iVertex] = latVertex[vertexIndices] - grid_corner_lon[cellIndices, iVertex] = lonVertex[vertexIndices] - - outFile.variables['grid_corner_lat'][:] = grid_corner_lat[:] - outFile.variables['grid_corner_lon'][:] = grid_corner_lon[:] - - # Update history attribute of netCDF file - if hasattr(inFile, 'history'): - newhist = '\n'.join([getattr(inFile, 'history'), - ' '.join(sys.argv[:])]) - else: - newhist = sys.argv[:] - setattr(outFile, 'history', newhist) - - inFile.close() - outFile.close() - - def to_esmf(self, esmfFileName): - """ - Create an ESMF mesh file for the mesh - - Parameters - ---------- - esmfFileName : str - The path to which the ESMF mesh file should be written - """ - with xarray.open_dataset(self.fileName) as ds: - - nodeCount = ds.sizes['nVertices'] - elementCount = ds.sizes['nCells'] - coordDim = 2 - - nodeCoords = numpy.zeros((nodeCount, coordDim)) - nodeCoords[:, 0] = ds.lonVertex.values - nodeCoords[:, 1] = ds.latVertex.values - centerCoords = numpy.zeros((elementCount, coordDim)) - centerCoords[:, 0] = ds.lonCell.values - centerCoords[:, 1] = ds.latCell.values - - elementConn = ds.verticesOnCell.values - elementConn[elementConn == nodeCount + 1] = -1 - - ds_out = xarray.Dataset() - ds_out['nodeCoords'] = (('nodeCount', 'coordDim'), nodeCoords) - ds_out.nodeCoords.attrs['units'] = 'radians' - ds_out['centerCoords'] = \ - (('elementCount', 'coordDim'), centerCoords) - ds_out.centerCoords.attrs['units'] = 'radians' - ds_out['elementConn'] = \ - (('elementCount', 'maxNodePElement'), elementConn) - ds_out.elementConn.attrs['start_index'] = 1 - ds_out.elementConn.attrs['_FillValue'] = -1 - ds_out['numElementConn'] = \ - (('elementCount',), ds.nEdgesOnCell.values) - ds_out.attrs['gridType'] = 'unstructured mesh' - ds_out.attrs['version'] = '0.9' + warnings.warn('MpasMeshDescriptor is deprecated and will be removed ' + 'in the next release. Use MpasCellMeshDescriptor ' + 'instead.', DeprecationWarning) - ds_out.to_netcdf(esmfFileName, format=self.format, - engine=self.engine) + super().__init__(fileName, meshName=meshName, vertices=vertices) From 820c766a66a57bec879f2741c39f32f57f376a67 Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Sat, 16 Sep 2023 22:18:25 +0200 Subject: [PATCH 2/6] Add MpasVertexMeshDescriptor --- pyremap/__init__.py | 2 + pyremap/descriptor/__init__.py | 3 + .../descriptor/mpas_vertex_mesh_descriptor.py | 159 ++++++++++++++++++ 3 files changed, 164 insertions(+) create mode 100644 pyremap/descriptor/mpas_vertex_mesh_descriptor.py diff --git a/pyremap/__init__.py b/pyremap/__init__.py index 2293536..e63bdd0 100644 --- a/pyremap/__init__.py +++ b/pyremap/__init__.py @@ -1,8 +1,10 @@ from pyremap.descriptor import ( LatLon2DGridDescriptor, LatLonGridDescriptor, + MpasCellMeshDescriptor, MpasEdgeMeshDescriptor, MpasMeshDescriptor, + MpasVertexMeshDescriptor, PointCollectionDescriptor, ProjectionGridDescriptor, get_lat_lon_descriptor, diff --git a/pyremap/descriptor/__init__.py b/pyremap/descriptor/__init__.py index df40e1b..a207146 100644 --- a/pyremap/descriptor/__init__.py +++ b/pyremap/descriptor/__init__.py @@ -20,6 +20,9 @@ from pyremap.descriptor.mpas_cell_mesh_descriptor import MpasCellMeshDescriptor from pyremap.descriptor.mpas_edge_mesh_descriptor import MpasEdgeMeshDescriptor from pyremap.descriptor.mpas_mesh_descriptor import MpasMeshDescriptor +from pyremap.descriptor.mpas_vertex_mesh_descriptor import ( + MpasVertexMeshDescriptor, +) from pyremap.descriptor.point_collection_descriptor import ( PointCollectionDescriptor, ) diff --git a/pyremap/descriptor/mpas_vertex_mesh_descriptor.py b/pyremap/descriptor/mpas_vertex_mesh_descriptor.py new file mode 100644 index 0000000..536f99d --- /dev/null +++ b/pyremap/descriptor/mpas_vertex_mesh_descriptor.py @@ -0,0 +1,159 @@ +# This software is open source software available under the BSD-3 license. +# +# Copyright (c) 2022 Triad National Security, LLC. All rights reserved. +# Copyright (c) 2022 Lawrence Livermore National Security, LLC. All rights +# reserved. +# Copyright (c) 2022 UT-Battelle, LLC. All rights reserved. +# +# Additional copyright and license information can be found in the LICENSE file +# distributed with this code, or at +# https://raw.githubusercontent.com/MPAS-Dev/pyremap/main/LICENSE + +import sys + +import netCDF4 +import numpy as np +import xarray as xr + +from pyremap.descriptor.mesh_descriptor import MeshDescriptor +from pyremap.descriptor.utility import create_scrip + + +class MpasVertexMeshDescriptor(MeshDescriptor): + """ + A class for describing an MPAS vertex mesh + + Attributes + ---------- + fileName : str + The path of the file containing the MPAS mesh + """ + def __init__(self, fileName, meshName=None): + """ + Constructor stores the file name + + Parameters + ---------- + fileName : str + The path of the file containing the MPAS mesh + + meshName : str, optional + The name of the MPAS mesh (e.g. ``'oEC60to30'`` or + ``'oRRS18to6'``). If not provided, the data set in ``fileName`` + must have a global attribute ``meshName`` that will be used + instead. + """ + super().__init__() + + with xr.open_dataset(fileName) as ds: + + if meshName is None: + if 'meshName' not in ds.attrs: + raise ValueError('No meshName provided or found in file.') + self.meshName = ds.attrs['meshName'] + else: + self.meshName = meshName + + self.fileName = fileName + self.regional = True + + # build coords + self.coords = {'latVertex': {'dims': 'nVertices', + 'data': ds.latVertex.values, + 'attrs': {'units': 'radians'}}, + 'lonVertex': {'dims': 'nVertices', + 'data': ds.lonVertex.values, + 'attrs': {'units': 'radians'}}} + self.dims = ['nVertices'] + self.dimSize = [ds.dims[dim] for dim in self.dims] + + def to_scrip(self, scripFileName): + """ + Given an MPAS mesh file, create a SCRIP file based on the mesh. + + Parameters + ---------- + scripFileName : str + The path to which the SCRIP file should be written + """ + + inFile = netCDF4.Dataset(self.fileName, 'r') + outFile = netCDF4.Dataset(scripFileName, 'w', format=self.format) + + # Get info from input file + latCell = inFile.variables['latCell'][:] + lonCell = inFile.variables['lonCell'][:] + latVertex = inFile.variables['latVertex'][:] + lonVertex = inFile.variables['lonVertex'][:] + latEdge = inFile.variables['latEdge'][:] + lonEdge = inFile.variables['lonEdge'][:] + cellsOnVertex = inFile.variables['cellsOnVertex'][:] + edgesOnVertex = inFile.variables['edgesOnVertex'][:] + nVertices = len(inFile.dimensions['nVertices']) + vertexDegree = len(inFile.dimensions['vertexDegree']) + kiteAreasOnVertex = inFile.variables['kiteAreasOnVertex'][:] + sphereRadius = float(inFile.sphere_radius) + + if vertexDegree != 3: + raise ValueError(f'MpasVertexMeshDescriptor does not support ' + f'vertexDegree {vertexDegree}') + + create_scrip(outFile, grid_size=nVertices, grid_corners=6, + grid_rank=1, units='radians', meshName=self.meshName) + + validCellsOnVertex = np.zeros(nVertices, dtype=int) + vertexArea = np.zeros(nVertices) + for iCell in range(vertexDegree): + mask = cellsOnVertex[:, iCell] > 0 + validCellsOnVertex[mask] = validCellsOnVertex[mask] + 1 + vertexArea[mask] = (vertexArea[mask] + + kiteAreasOnVertex[mask, iCell]) + + grid_area = outFile.createVariable('grid_area', 'f8', ('grid_size',)) + grid_area.units = 'radian^2' + # SCRIP uses square radians + grid_area[:] = vertexArea[:] / (sphereRadius**2) + + outFile.variables['grid_center_lat'][:] = latVertex[:] + outFile.variables['grid_center_lon'][:] = lonVertex[:] + outFile.variables['grid_dims'][:] = nVertices + outFile.variables['grid_imask'][:] = 1 + + # grid corners: + grid_corner_lon = np.zeros((nVertices, 6)) + grid_corner_lat = np.zeros((nVertices, 6)) + + # start by filling the corners with the vertex lat and lon as the + # fallback + for iCorner in range(6): + grid_corner_lon[:, iCorner] = lonVertex + grid_corner_lat[:, iCorner] = latVertex + + # if edges on vertex exist, replace even indices with their locations + for iEdge in range(3): + mask = edgesOnVertex[:, iEdge] > 0 + edges = edgesOnVertex[mask, iEdge] - 1 + grid_corner_lon[mask, 2 * iEdge] = lonEdge[edges] + grid_corner_lat[mask, 2 * iEdge] = latEdge[edges] + + # similarly, if cells on vertex exist, replace odd indices with their + # locations + for iCell in range(3): + mask = cellsOnVertex[:, iCell] > 0 + cells = cellsOnVertex[mask, iCell] - 1 + grid_corner_lon[mask, 2 * iCell + 1] = lonCell[cells] + grid_corner_lat[mask, 2 * iCell + 1] = latCell[cells] + + outFile.variables['grid_corner_lat'][:] = grid_corner_lat[:] + outFile.variables['grid_corner_lon'][:] = grid_corner_lon[:] + + # Update history attribute of netCDF file + if hasattr(inFile, 'history'): + newhist = '\n'.join([getattr(inFile, 'history'), + ' '.join(sys.argv[:])]) + else: + newhist = sys.argv[:] + setattr(outFile, 'history', newhist) + + inFile.close() + outFile.close() From c999a700fdc4230a36c002db31f139bce526acea Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Mon, 18 Sep 2023 08:42:11 +0200 Subject: [PATCH 3/6] Update remapper for new descriptors --- pyremap/remapper.py | 199 ++++++++++++++++++-------------------------- 1 file changed, 80 insertions(+), 119 deletions(-) diff --git a/pyremap/remapper.py b/pyremap/remapper.py index b50430f..ac02d53 100644 --- a/pyremap/remapper.py +++ b/pyremap/remapper.py @@ -25,8 +25,9 @@ from pyremap.descriptor import ( LatLon2DGridDescriptor, LatLonGridDescriptor, + MpasCellMeshDescriptor, MpasEdgeMeshDescriptor, - MpasMeshDescriptor, + MpasVertexMeshDescriptor, PointCollectionDescriptor, ProjectionGridDescriptor, ) @@ -74,14 +75,14 @@ def __init__(self, sourceDescriptor, destinationDescriptor, raise TypeError("sourceDescriptor of type " "PointCollectionDescriptor is not supported.") if not isinstance(sourceDescriptor, - (MpasMeshDescriptor, MpasEdgeMeshDescriptor, - LatLonGridDescriptor, + (MpasCellMeshDescriptor, MpasEdgeMeshDescriptor, + MpasVertexMeshDescriptor, LatLonGridDescriptor, LatLon2DGridDescriptor, ProjectionGridDescriptor)): raise TypeError("sourceDescriptor is not of a recognized type.") if not isinstance(destinationDescriptor, - (MpasMeshDescriptor, MpasEdgeMeshDescriptor, - LatLonGridDescriptor, + (MpasCellMeshDescriptor, MpasEdgeMeshDescriptor, + MpasVertexMeshDescriptor, LatLonGridDescriptor, LatLon2DGridDescriptor, ProjectionGridDescriptor, PointCollectionDescriptor)): raise TypeError( @@ -96,7 +97,8 @@ def __init__(self, sourceDescriptor, destinationDescriptor, def build_mapping_file(self, method='bilinear', # noqa: C901 additionalArgs=None, logger=None, mpiTasks=1, tempdir=None, esmf_path=None, - esmf_parallel_exec=None, extrap_method=None): + esmf_parallel_exec=None, extrap_method=None, + include_logs=False): """ Given a source file defining either an MPAS mesh or a lat-lon grid and a destination file or set of arrays defining a lat-lon grid, constructs @@ -137,6 +139,9 @@ def build_mapping_file(self, method='bilinear', # noqa: C901 extrap_method : {'neareststod', 'nearestidavg','creep'}, optional The method used to extrapolate unmapped destination locations + include_logs : bool, optional + Whether to include log files from ``ESMF_RegridWeightGen`` + Raises ------ OSError @@ -152,9 +157,8 @@ def build_mapping_file(self, method='bilinear', # noqa: C901 if isinstance(self.destinationDescriptor, PointCollectionDescriptor) and \ method not in ['bilinear', 'neareststod']: - raise ValueError("method {} not supported for destination " - "grid of type PointCollectionDescriptor." - "".format(method)) + raise ValueError(f'method {method} not supported for destination ' + f'grid of type PointCollectionDescriptor.') if self.mappingFileName is None or \ os.path.exists(self.mappingFileName): @@ -182,67 +186,26 @@ def build_mapping_file(self, method='bilinear', # noqa: C901 else: tempobj = None - sourceFileName = '{}/src_mesh.nc'.format(tempdir) - destinationFileName = '{}/dst_mesh.nc'.format(tempdir) + sourceFileName = f'{tempdir}/src_mesh.nc' + destinationFileName = f'{tempdir}/dst_mesh.nc' - src_loc = 'center' - src_file_format = 'scrip' - if isinstance(self.sourceDescriptor, - (MpasMeshDescriptor, MpasEdgeMeshDescriptor)): - src_file_format = 'esmf' + self.sourceDescriptor.to_scrip(sourceFileName) - if isinstance(self.sourceDescriptor, MpasMeshDescriptor) and \ - self.sourceDescriptor.vertices: - if 'conserve' in method: - raise ValueError('Can\'t remap from MPAS vertices with ' - 'conservative methods') - src_loc = 'corner' - - dst_loc = 'center' - dst_file_format = 'scrip' - if isinstance(self.destinationDescriptor, - (MpasMeshDescriptor, MpasEdgeMeshDescriptor)): - dst_file_format = 'esmf' - - if isinstance(self.destinationDescriptor, MpasMeshDescriptor) and \ - self.destinationDescriptor.vertices: - if 'conserve' in method: - raise ValueError('Can\'t remap to MPAS vertices with ' - 'conservative methods') - dst_loc = 'corner' - - if src_file_format == 'scrip': - self.sourceDescriptor.to_scrip(sourceFileName) - elif src_file_format == 'esmf': - self.sourceDescriptor.to_esmf(sourceFileName) - else: - raise ValueError('Unexpected file format {}'.format( - src_file_format)) - - if dst_file_format == 'scrip': - self.destinationDescriptor.to_scrip(destinationFileName) - elif dst_file_format == 'esmf': - self.destinationDescriptor.to_esmf(destinationFileName) - else: - raise ValueError('Unexpected file format {}'.format( - dst_file_format)) + self.destinationDescriptor.to_scrip(destinationFileName) args = [rwgPath, '--source', sourceFileName, '--destination', destinationFileName, '--weight', self.mappingFileName, '--method', method, - '--netcdf4', - '--no_log'] + '--netcdf4'] + + if not include_logs: + args.extend(['--no_log']) if extrap_method is not None: args.extend(['--extrap_method', extrap_method]) - if src_file_format == 'esmf': - args.extend(['--src_loc', src_loc]) - if dst_file_format == 'esmf': - args.extend(['--dst_loc', dst_loc]) - parallel_args = [] if esmf_parallel_exec is not None: @@ -250,10 +213,10 @@ def build_mapping_file(self, method='bilinear', # noqa: C901 parallel_args = esmf_parallel_exec.split(' ') if 'srun' in esmf_parallel_exec: - parallel_args.extend(['-n', '{}'.format(mpiTasks)]) + parallel_args.extend(['-n', f'{mpiTasks}']) else: # presume mpirun syntax - parallel_args.extend(['-np', '{}'.format(mpiTasks)]) + parallel_args.extend(['-np', f'{mpiTasks}']) elif 'CONDA_PREFIX' in os.environ and mpiTasks > 1: # this is a conda environment, so we need to find out if esmf @@ -265,14 +228,13 @@ def build_mapping_file(self, method='bilinear', # noqa: C901 if 'mpi_mpich' in build_string or 'mpi_openmpi' in build_string: # esmf was installed with MPI, so we should use mpirun - mpirun_path = '{}/bin/mpirun'.format( - os.environ['CONDA_PREFIX']) - parallel_args = [mpirun_path, '-np', '{}'.format(mpiTasks)] + mpirun_path = f'{os.environ["CONDA_PREFIX"]}/bin/mpirun' + parallel_args = [mpirun_path, '-np', f'{mpiTasks}'] else: # esmf was installed without MPI, so we shouldn't try to # use it - warnings.warn('Requesting {} MPI tasks but the MPI version' - ' of ESMF is not installed'.format(mpiTasks)) + warnings.warn(f'Requesting {mpiTasks} MPI tasks but the MPI ' + f'version of ESMF is not installed') args = parallel_args + args @@ -379,10 +341,10 @@ def remap_file(self, inFileName, outFileName, # noqa: C901 # Xylar Asay-Davis if self.mappingFileName is None: - raise ValueError('No mapping file was given because remapping is ' - 'not necessary. The calling\n' - 'code should simply use the constents of {} ' - 'directly.'.format(inFileName)) + raise ValueError(f'No mapping file was given because remapping is ' + f'not necessary. The calling\n' + f'code should simply use the constents of ' + f'{inFileName} directly.') if not overwrite and os.path.exists(outFileName): # a remapped file already exists, so nothing to do @@ -411,22 +373,22 @@ def remap_file(self, inFileName, outFileName, # noqa: C901 regridArgs = [] - if isinstance(self.sourceDescriptor, MpasMeshDescriptor): - if self.sourceDescriptor.vertices: - regridArgs.extend(['--rgr col_nm=nVertices']) - else: - args.extend(['-P', 'mpas']) - if not replaceMpasFill: - # the -C (climatology) flag prevents ncremap from trying to - # add a _FillValue attribute that might already be present - # and quits with an error - args.append('-C') + if isinstance(self.sourceDescriptor, (MpasCellMeshDescriptor, + MpasEdgeMeshDescriptor, + MpasVertexMeshDescriptor)): + args.extend(['-P', 'mpas']) + if not replaceMpasFill: + # the -C (climatology) flag prevents ncremap from trying to + # add a _FillValue attribute that might already be present + # and quits with an error + args.append('-C') if isinstance(self.sourceDescriptor, MpasEdgeMeshDescriptor): regridArgs.extend(['--rgr col_nm=nEdges']) - if isinstance(self.sourceDescriptor, - (MpasMeshDescriptor, MpasEdgeMeshDescriptor)) and \ + if isinstance(self.sourceDescriptor, (MpasCellMeshDescriptor, + MpasEdgeMeshDescriptor, + MpasVertexMeshDescriptor)) and \ renormalize is not None: # we also want to make sure cells that receive no data are # marked with fill values, even if the source MPAS data @@ -437,31 +399,28 @@ def remap_file(self, inFileName, outFileName, # noqa: C901 args.extend(['-v', ','.join(variableList)]) if renormalize is not None: - regridArgs.append('--renormalize={}'.format(renormalize)) + regridArgs.append(f'--renormalize={renormalize}') if isinstance(self.sourceDescriptor, LatLonGridDescriptor): - regridArgs.extend(['--rgr lat_nm={}'.format( - self.sourceDescriptor.latVarName), - '--rgr lon_nm={}'.format( - self.sourceDescriptor.lonVarName)]) + regridArgs.extend( + [f'--rgr lat_nm={self.sourceDescriptor.latVarName}', + f'--rgr lon_nm={self.sourceDescriptor.lonVarName}']) elif isinstance(self.sourceDescriptor, ProjectionGridDescriptor): - regridArgs.extend(['--rgr lat_nm={}'.format( - self.sourceDescriptor.yVarName), - '--rgr lon_nm={}'.format( - self.sourceDescriptor.xVarName)]) + regridArgs.extend( + [f'--rgr lat_nm={self.sourceDescriptor.yVarName}', + f'--rgr lon_nm={self.sourceDescriptor.xVarName}']) if isinstance(self.destinationDescriptor, LatLonGridDescriptor): - regridArgs.extend([ - f'--rgr lat_nm_out={self.destinationDescriptor.latVarName}', - f'--rgr lon_nm_out={self.destinationDescriptor.lonVarName}', - f'--rgr lat_dmn_nm={self.destinationDescriptor.latVarName}', - f'--rgr lon_dmn_nm={self.destinationDescriptor.lonVarName}']) + regridArgs.extend( + [f'--rgr lat_nm_out={self.destinationDescriptor.latVarName}', + f'--rgr lon_nm_out={self.destinationDescriptor.lonVarName}', + f'--rgr lat_dmn_nm={self.destinationDescriptor.latVarName}', + f'--rgr lon_dmn_nm={self.destinationDescriptor.lonVarName}']) elif isinstance(self.destinationDescriptor, ProjectionGridDescriptor): - regridArgs.extend(['--rgr lat_dmn_nm={}'.format( - self.destinationDescriptor.yVarName), - '--rgr lon_dmn_nm={}'.format( - self.destinationDescriptor.xVarName), - '--rgr lat_nm_out=lat', '--rgr lon_nm_out=lon']) + regridArgs.extend( + [f'--rgr lat_dmn_nm={self.destinationDescriptor.yVarName}', + f'--rgr lon_dmn_nm={self.destinationDescriptor.xVarName}', + '--rgr lat_nm_out=lat', '--rgr lon_nm_out=lon']) if isinstance(self.destinationDescriptor, PointCollectionDescriptor): regridArgs.extend(['--rgr lat_nm_out=lat', '--rgr lon_nm_out=lon']) @@ -550,10 +509,10 @@ def remap(self, ds, renormalizationThreshold=None): for index, dim in enumerate(self.sourceDescriptor.dims): if self.src_grid_dims[index] != ds.sizes[dim]: - raise ValueError('data set and remapping source dimension {} ' - 'don\'t have the same size: {} != {}'.format( - dim, self.src_grid_dims[index], - ds.sizes[dim])) + raise ValueError( + f'data set and remapping source dimension {dim} don\'t ' + f'have the same size: {self.src_grid_dims[index]} != ' + f'{ds.sizes[dim]}') if isinstance(ds, xr.DataArray): remappedDs = self._remap_data_array(ds, renormalizationThreshold) @@ -605,13 +564,13 @@ def _load_mapping(self): # check that the mapping file has the right number of dimensions if nSourceDims != src_grid_rank or \ nDestinationDims != dst_grid_rank: - raise ValueError('The number of source and/or ' - 'destination dimensions does not\n' - 'match the expected number of source and ' - 'destination dimensions in the mapping\n' - 'file. {} != {} and/or {} != {}'.format( - nSourceDims, src_grid_rank, - nDestinationDims, dst_grid_rank)) + raise ValueError( + f'The number of source and/or destination dimensions does not ' + f'match the expected \n' + f'number of source and destination dimensions in the mapping ' + f'file. \n' + f'{nSourceDims} != {src_grid_rank} and/or {nDestinationDims} ' + f'!= {dst_grid_rank}') # grid dimensions need to be reversed because they are in Fortran order self.src_grid_dims = dsMapping['src_grid_dims'].values[::-1] @@ -623,17 +582,19 @@ def _load_mapping(self): dimSize = self.sourceDescriptor.dimSize[index] checkDimSize = self.src_grid_dims[index] if dimSize != checkDimSize: - raise ValueError('source mesh descriptor and remapping source ' - 'dimension {} don\'t have the same size: \n' - '{} != {}'.format(dim, dimSize, checkDimSize)) + raise ValueError( + f'source mesh descriptor and remapping source dimension ' + f'{dim} don\'t have the same size: \n' + f'{dimSize} != {checkDimSize}') for index in range(len(self.destinationDescriptor.dims)): dim = self.destinationDescriptor.dims[index] dimSize = self.destinationDescriptor.dimSize[index] checkDimSize = self.dst_grid_dims[index] if dimSize != checkDimSize: - raise ValueError('dest. mesh descriptor and remapping dest. ' - 'dimension {} don\'t have the same size: \n' - '{} != {}'.format(dim, dimSize, checkDimSize)) + raise ValueError( + f'dest. mesh descriptor and remapping dest. dimension ' + f'{dim} don\'t have the same size: \n' + f'{dimSize} != {checkDimSize}') self.frac_b = dsMapping['frac_b'].values @@ -791,6 +752,6 @@ def _print_running(args, fn): print_args = [] for arg in args: if ' ' in arg: - arg = '"{}"'.format(arg) + arg = f'"{arg}"' print_args.append(arg) - fn('running: {}'.format(' '.join(print_args))) + fn(f'running: {" ".join(print_args)}') From 8b3640cc424dfd677cdb0c614c3b50f3c9c85767 Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Mon, 18 Sep 2023 08:24:31 +0200 Subject: [PATCH 4/6] Update examples --- .../make_mpas_to_Antarctic_stereo_mapping.py | 16 ++--- examples/make_mpas_to_lat_lon_mapping.py | 12 ++-- ...mpas_vertex_to_Antarctic_stereo_mapping.py | 66 +++++++++++++++++++ examples/remap_stereographic.py | 13 ++-- 4 files changed, 86 insertions(+), 21 deletions(-) create mode 100755 examples/make_mpas_vertex_to_Antarctic_stereo_mapping.py diff --git a/examples/make_mpas_to_Antarctic_stereo_mapping.py b/examples/make_mpas_to_Antarctic_stereo_mapping.py index 681aeed..6eafa12 100755 --- a/examples/make_mpas_to_Antarctic_stereo_mapping.py +++ b/examples/make_mpas_to_Antarctic_stereo_mapping.py @@ -10,18 +10,18 @@ # distributed with this code, or at # https://raw.githubusercontent.com/MPAS-Dev/pyremap/main/LICENSE -''' +""" Creates a mapping file that can be used with ncremap (NCO) to remap MPAS files to a latitude/longitude grid. Usage: Copy this script into the main MPAS-Analysis directory (up one level). Modify the grid name, the path to the MPAS grid file and the output grid resolution. -''' +""" import xarray -from pyremap import MpasMeshDescriptor, Remapper, get_polar_descriptor +from pyremap import MpasCellMeshDescriptor, Remapper, get_polar_descriptor # replace with the MPAS mesh name @@ -32,14 +32,14 @@ # https://web.lcrc.anl.gov/public/e3sm/inputdata/ocn/mpas-o/oQU240/ocean.QU.240km.151209.nc inGridFileName = 'ocean.QU.240km.151209.nc' -inDescriptor = MpasMeshDescriptor(inGridFileName, inGridName) +inDescriptor = MpasCellMeshDescriptor(inGridFileName, inGridName) # modify the size and resolution of the Antarctic grid as desired outDescriptor = get_polar_descriptor(Lx=6000., Ly=5000., dx=10., dy=10., projection='antarctic') outGridName = outDescriptor.meshName -mappingFileName = 'map_{}_to_{}_bilinear.nc'.format(inGridName, outGridName) +mappingFileName = f'map_{inGridName}_to_{outGridName}_bilinear.nc' remapper = Remapper(inDescriptor, outDescriptor, mappingFileName) @@ -47,17 +47,17 @@ remapper.build_mapping_file(method='bilinear', mpiTasks=4) # select the SST at the initial time as an example data set -srcFileName = 'temp_{}.nc'.format(inGridName) +srcFileName = f'temp_{inGridName}.nc' ds = xarray.open_dataset(inGridFileName) dsOut = xarray.Dataset() dsOut['temperature'] = ds['temperature'].isel(nVertLevels=0, Time=0) dsOut.to_netcdf(srcFileName) # do remapping with ncremap -outFileName = 'temp_{}_file.nc'.format(outGridName) +outFileName = f'temp_{outGridName}_file.nc' remapper.remap_file(srcFileName, outFileName) # do remapping again, this time with python remapping -outFileName = 'temp_{}_array.nc'.format(outGridName) +outFileName = f'temp_{outGridName}_array.nc' dsOut = remapper.remap(dsOut) dsOut.to_netcdf(outFileName) diff --git a/examples/make_mpas_to_lat_lon_mapping.py b/examples/make_mpas_to_lat_lon_mapping.py index ead373b..c12784f 100755 --- a/examples/make_mpas_to_lat_lon_mapping.py +++ b/examples/make_mpas_to_lat_lon_mapping.py @@ -10,18 +10,18 @@ # distributed with this code, or at # https://raw.githubusercontent.com/MPAS-Dev/pyremap/main/LICENSE -''' +""" Creates a mapping file that can be used with ncremap (NCO) to remap MPAS files to a latitude/longitude grid. Usage: Copy this script into the main MPAS-Analysis directory (up one level). Modify the grid name, the path to the MPAS grid file and the output grid resolution. -''' +""" import xarray -from pyremap import MpasMeshDescriptor, Remapper, get_lat_lon_descriptor +from pyremap import MpasCellMeshDescriptor, Remapper, get_lat_lon_descriptor # replace with the MPAS mesh name inGridName = 'oQU240' @@ -31,20 +31,20 @@ # https://web.lcrc.anl.gov/public/e3sm/inputdata/ocn/mpas-o/oQU240/ocean.QU.240km.151209.nc inGridFileName = 'ocean.QU.240km.151209.nc' -inDescriptor = MpasMeshDescriptor(inGridFileName, inGridName) +inDescriptor = MpasCellMeshDescriptor(inGridFileName, inGridName) # modify the resolution of the global lat-lon grid as desired outDescriptor = get_lat_lon_descriptor(dLon=0.5, dLat=0.5) outGridName = outDescriptor.meshName -mappingFileName = 'map_{}_to_{}_bilinear.nc'.format(inGridName, outGridName) +mappingFileName = f'map_{inGridName}_to_{outGridName}_bilinear.nc' remapper = Remapper(inDescriptor, outDescriptor, mappingFileName) remapper.build_mapping_file(method='bilinear', mpiTasks=1) -outFileName = 'temp_{}.nc'.format(outGridName) +outFileName = f'temp_{outGridName}.nc' ds = xarray.open_dataset(inGridFileName) dsOut = xarray.Dataset() dsOut['temperature'] = ds['temperature'] diff --git a/examples/make_mpas_vertex_to_Antarctic_stereo_mapping.py b/examples/make_mpas_vertex_to_Antarctic_stereo_mapping.py new file mode 100755 index 0000000..cfa81ad --- /dev/null +++ b/examples/make_mpas_vertex_to_Antarctic_stereo_mapping.py @@ -0,0 +1,66 @@ +#!/usr/bin/env python +# This software is open source software available under the BSD-3 license. +# +# Copyright (c) 2019 Triad National Security, LLC. All rights reserved. +# Copyright (c) 2019 Lawrence Livermore National Security, LLC. All rights +# reserved. +# Copyright (c) 2019 UT-Battelle, LLC. All rights reserved. +# +# Additional copyright and license information can be found in the LICENSE file +# distributed with this code, or at +# https://raw.githubusercontent.com/MPAS-Dev/pyremap/main/LICENSE + +""" +Creates a mapping file that can be used with ncremap (NCO) to remap MPAS files +to a latitude/longitude grid. + +Usage: Copy this script into the main MPAS-Analysis directory (up one level). +Modify the grid name, the path to the MPAS grid file and the output grid +resolution. +""" + +import numpy as np +import xarray as xr + +from pyremap import MpasVertexMeshDescriptor, Remapper, get_polar_descriptor + + +# replace with the MPAS mesh name +inGridName = 'oQU240' + +# replace with the path to the desired mesh or restart file +# As an example, use: +# https://web.lcrc.anl.gov/public/e3sm/inputdata/ocn/mpas-o/oQU240/ocean.QU.240km.151209.nc +inGridFileName = 'ocean.QU.240km.151209.nc' + +inDescriptor = MpasVertexMeshDescriptor(inGridFileName, inGridName) + +# modify the size and resolution of the Antarctic grid as desired +outDescriptor = get_polar_descriptor(Lx=6000., Ly=5000., dx=10., dy=10., + projection='antarctic') +outGridName = outDescriptor.meshName + +mappingFileName = f'map_{inGridName}_vertex_to_{outGridName}_conserve.nc' + +remapper = Remapper(inDescriptor, outDescriptor, mappingFileName) + +# conservative remapping with 4 MPI tasks (using mpirun) +remapper.build_mapping_file(method='conserve', mpiTasks=4, + esmf_parallel_exec='mpirun', include_logs=True) + +# select the SST at the initial time as an example data set +srcFileName = f'vertex_id_{inGridName}.nc' +ds = xr.open_dataset(inGridFileName) +dsOut = xr.Dataset() +dsOut['indexToVertexID'] = ds['indexToVertexID'].astype(float) +dsOut['random'] = ('nVertices', np.random.random(ds.sizes['nVertices'])) +dsOut.to_netcdf(srcFileName) + +# do remapping with ncremap +outFileName = f'vertex_id_{outGridName}_file.nc' +remapper.remap_file(srcFileName, outFileName) + +# do remapping again, this time with python remapping +outFileName = f'vertex_id_{outGridName}_array.nc' +dsOut = remapper.remap(dsOut) +dsOut.to_netcdf(outFileName) diff --git a/examples/remap_stereographic.py b/examples/remap_stereographic.py index ff5cd6c..2e1af13 100755 --- a/examples/remap_stereographic.py +++ b/examples/remap_stereographic.py @@ -1,13 +1,13 @@ #!/usr/bin/env python -''' +""" Creates a mapping file and remaps the variables from an input file on an Antarctic grid to another grid with the same extent but a different resolution. The mapping file can be used with ncremap (NCO) to remap MPAS files between these same gridss. Usage: Copy this script into the main MPAS-Analysis directory (up one level). -''' +""" import numpy import xarray @@ -33,7 +33,7 @@ if args.method not in ['bilinear', 'neareststod', 'conserve']: - raise ValueError('Unexpected method {}'.format(args.method)) + raise ValueError(f'Unexpected method {args.method}') dsIn = xarray.open_dataset(args.inFileName) @@ -43,7 +43,7 @@ Lx = int((x[-1] - x[0])/1000.) Ly = int((y[-1] - y[0])/1000.) -inMeshName = '{}x{}km_{}km_Antarctic_stereo'.format(Lx, Ly, dx) +inMeshName = f'{Lx}x{Ly}km_{dx}km_Antarctic_stereo' projection = get_antarctic_stereographic_projection() @@ -58,13 +58,12 @@ yOut = y[0] + outRes*numpy.arange(nyOut) -outMeshName = '{}x{}km_{}km_Antarctic_stereo'.format(Lx, Ly, args.resolution) +outMeshName = f'{Lx}x{Ly}km_{args.resolution}km_Antarctic_stereo' outDescriptor = ProjectionGridDescriptor.create(projection, xOut, yOut, outMeshName) -mappingFileName = 'map_{}_to_{}_{}.nc'.format(inMeshName, outMeshName, - args.method) +mappingFileName = f'map_{inMeshName}_to_{outMeshName}_{args.method}.nc' remapper = Remapper(inDescriptor, outDescriptor, mappingFileName) From 4ee818366c6eec01727d1063ccde03b0963e2674 Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Mon, 18 Sep 2023 12:09:07 +0200 Subject: [PATCH 5/6] Write history as string, not string array --- .../descriptor/lat_lon_2d_grid_descriptor.py | 9 ++------- pyremap/descriptor/lat_lon_grid_descriptor.py | 11 +++-------- pyremap/descriptor/mpas_cell_mesh_descriptor.py | 16 +++++++--------- pyremap/descriptor/mpas_edge_mesh_descriptor.py | 17 +++++++---------- .../descriptor/mpas_vertex_mesh_descriptor.py | 17 +++++++---------- .../descriptor/point_collection_descriptor.py | 10 ++++++---- .../descriptor/projection_grid_descriptor.py | 16 +++++++--------- pyremap/descriptor/utility.py | 15 ++++++++++++++- 8 files changed, 53 insertions(+), 58 deletions(-) diff --git a/pyremap/descriptor/lat_lon_2d_grid_descriptor.py b/pyremap/descriptor/lat_lon_2d_grid_descriptor.py index 0c38596..8e8727a 100644 --- a/pyremap/descriptor/lat_lon_2d_grid_descriptor.py +++ b/pyremap/descriptor/lat_lon_2d_grid_descriptor.py @@ -9,14 +9,13 @@ # distributed with this code, or at # https://raw.githubusercontent.com/MPAS-Dev/pyremap/main/LICENSE -import sys - import netCDF4 import numpy import xarray from pyremap.descriptor.mesh_descriptor import MeshDescriptor from pyremap.descriptor.utility import ( + add_history, create_scrip, interp_extrap_corners_2d, round_res, @@ -120,11 +119,7 @@ def read(cls, fileName=None, ds=None, latVarName='lat', descriptor._set_coords(latVarName, lonVarName, ds[latVarName].dims[0], ds[latVarName].dims[1]) - if 'history' in ds.attrs: - descriptor.history = '\n'.join([ds.attrs['history'], - ' '.join(sys.argv[:])]) - else: - descriptor.history = sys.argv[:] + descriptor.history = add_history(ds=ds) return descriptor def to_scrip(self, scripFileName): diff --git a/pyremap/descriptor/lat_lon_grid_descriptor.py b/pyremap/descriptor/lat_lon_grid_descriptor.py index cc56ebb..6fb4d2b 100644 --- a/pyremap/descriptor/lat_lon_grid_descriptor.py +++ b/pyremap/descriptor/lat_lon_grid_descriptor.py @@ -9,14 +9,13 @@ # distributed with this code, or at # https://raw.githubusercontent.com/MPAS-Dev/pyremap/main/LICENSE -import sys - import netCDF4 import numpy import xarray from pyremap.descriptor.mesh_descriptor import MeshDescriptor from pyremap.descriptor.utility import ( + add_history, create_scrip, interp_extrap_corner, round_res, @@ -157,11 +156,7 @@ def read(cls, fileName=None, ds=None, latVarName='lat', descriptor._set_coords(latVarName, lonVarName, ds[latVarName].dims[0], ds[lonVarName].dims[0]) - if 'history' in ds.attrs: - descriptor.history = '\n'.join([ds.attrs['history'], - ' '.join(sys.argv[:])]) - else: - descriptor.history = sys.argv[:] + descriptor.history = add_history(ds=ds) return descriptor @classmethod @@ -199,7 +194,7 @@ def create(cls, latCorner, lonCorner, units='degrees', meshName=None, descriptor.lon = 0.5 * (lonCorner[0:-1] + lonCorner[1:]) descriptor.lat = 0.5 * (latCorner[0:-1] + latCorner[1:]) descriptor.units = units - descriptor.history = sys.argv[:] + descriptor.history = add_history() descriptor._set_coords('lat', 'lon', 'lat', 'lon') return descriptor diff --git a/pyremap/descriptor/mpas_cell_mesh_descriptor.py b/pyremap/descriptor/mpas_cell_mesh_descriptor.py index 71fac7b..904f76d 100644 --- a/pyremap/descriptor/mpas_cell_mesh_descriptor.py +++ b/pyremap/descriptor/mpas_cell_mesh_descriptor.py @@ -9,7 +9,6 @@ # distributed with this code, or at # https://raw.githubusercontent.com/MPAS-Dev/pyremap/main/LICENSE -import sys import warnings import netCDF4 @@ -17,7 +16,7 @@ import xarray from pyremap.descriptor.mesh_descriptor import MeshDescriptor -from pyremap.descriptor.utility import create_scrip +from pyremap.descriptor.utility import add_history, create_scrip class MpasCellMeshDescriptor(MeshDescriptor): @@ -32,6 +31,9 @@ class MpasCellMeshDescriptor(MeshDescriptor): fileName : str The path of the file containing the MPAS mesh + + history : str + The history attribute written to SCRIP files """ def __init__(self, fileName, meshName=None, vertices=False): """ @@ -94,6 +96,8 @@ def __init__(self, fileName, meshName=None, vertices=False): self.dims = ['nCells'] self.dimSize = [ds.dims[dim] for dim in self.dims] + self.history = add_history(ds=ds) + def to_scrip(self, scripFileName): """ Given an MPAS mesh file, create a SCRIP file based on the mesh. @@ -148,13 +152,7 @@ def to_scrip(self, scripFileName): outFile.variables['grid_corner_lat'][:] = grid_corner_lat[:] outFile.variables['grid_corner_lon'][:] = grid_corner_lon[:] - # Update history attribute of netCDF file - if hasattr(inFile, 'history'): - newhist = '\n'.join([getattr(inFile, 'history'), - ' '.join(sys.argv[:])]) - else: - newhist = sys.argv[:] - setattr(outFile, 'history', newhist) + setattr(outFile, 'history', self.history) inFile.close() outFile.close() diff --git a/pyremap/descriptor/mpas_edge_mesh_descriptor.py b/pyremap/descriptor/mpas_edge_mesh_descriptor.py index d2cbf3a..0d385c6 100644 --- a/pyremap/descriptor/mpas_edge_mesh_descriptor.py +++ b/pyremap/descriptor/mpas_edge_mesh_descriptor.py @@ -9,14 +9,12 @@ # distributed with this code, or at # https://raw.githubusercontent.com/MPAS-Dev/pyremap/main/LICENSE -import sys - import netCDF4 import numpy as np import xarray as xr from pyremap.descriptor.mesh_descriptor import MeshDescriptor -from pyremap.descriptor.utility import create_scrip +from pyremap.descriptor.utility import add_history, create_scrip class MpasEdgeMeshDescriptor(MeshDescriptor): @@ -27,6 +25,9 @@ class MpasEdgeMeshDescriptor(MeshDescriptor): ---------- fileName : str The path of the file containing the MPAS mesh + + history : str + The history attribute written to SCRIP files """ def __init__(self, fileName, meshName=None): """ @@ -67,6 +68,8 @@ def __init__(self, fileName, meshName=None): self.dims = ['nEdges'] self.dimSize = [ds.dims[dim] for dim in self.dims] + self.history = add_history(ds=ds) + def to_scrip(self, scripFileName): """ Given an MPAS mesh file, create a SCRIP file based on the mesh. @@ -143,13 +146,7 @@ def to_scrip(self, scripFileName): outFile.variables['grid_corner_lat'][:] = grid_corner_lat[:] outFile.variables['grid_corner_lon'][:] = grid_corner_lon[:] - # Update history attribute of netCDF file - if hasattr(inFile, 'history'): - newhist = '\n'.join([getattr(inFile, 'history'), - ' '.join(sys.argv[:])]) - else: - newhist = sys.argv[:] - setattr(outFile, 'history', newhist) + setattr(outFile, 'history', self.history) inFile.close() outFile.close() diff --git a/pyremap/descriptor/mpas_vertex_mesh_descriptor.py b/pyremap/descriptor/mpas_vertex_mesh_descriptor.py index 536f99d..0e6682f 100644 --- a/pyremap/descriptor/mpas_vertex_mesh_descriptor.py +++ b/pyremap/descriptor/mpas_vertex_mesh_descriptor.py @@ -9,14 +9,12 @@ # distributed with this code, or at # https://raw.githubusercontent.com/MPAS-Dev/pyremap/main/LICENSE -import sys - import netCDF4 import numpy as np import xarray as xr from pyremap.descriptor.mesh_descriptor import MeshDescriptor -from pyremap.descriptor.utility import create_scrip +from pyremap.descriptor.utility import add_history, create_scrip class MpasVertexMeshDescriptor(MeshDescriptor): @@ -27,6 +25,9 @@ class MpasVertexMeshDescriptor(MeshDescriptor): ---------- fileName : str The path of the file containing the MPAS mesh + + history : str + The history attribute written to SCRIP files """ def __init__(self, fileName, meshName=None): """ @@ -67,6 +68,8 @@ def __init__(self, fileName, meshName=None): self.dims = ['nVertices'] self.dimSize = [ds.dims[dim] for dim in self.dims] + self.history = add_history(ds=ds) + def to_scrip(self, scripFileName): """ Given an MPAS mesh file, create a SCRIP file based on the mesh. @@ -147,13 +150,7 @@ def to_scrip(self, scripFileName): outFile.variables['grid_corner_lat'][:] = grid_corner_lat[:] outFile.variables['grid_corner_lon'][:] = grid_corner_lon[:] - # Update history attribute of netCDF file - if hasattr(inFile, 'history'): - newhist = '\n'.join([getattr(inFile, 'history'), - ' '.join(sys.argv[:])]) - else: - newhist = sys.argv[:] - setattr(outFile, 'history', newhist) + setattr(outFile, 'history', self.history) inFile.close() outFile.close() diff --git a/pyremap/descriptor/point_collection_descriptor.py b/pyremap/descriptor/point_collection_descriptor.py index 2e8e248..679e7e3 100644 --- a/pyremap/descriptor/point_collection_descriptor.py +++ b/pyremap/descriptor/point_collection_descriptor.py @@ -9,14 +9,12 @@ # distributed with this code, or at # https://raw.githubusercontent.com/MPAS-Dev/pyremap/main/LICENSE -import sys - import netCDF4 import numpy import xarray from pyremap.descriptor.mesh_descriptor import MeshDescriptor -from pyremap.descriptor.utility import create_scrip +from pyremap.descriptor.utility import add_history, create_scrip class PointCollectionDescriptor(MeshDescriptor): @@ -33,6 +31,9 @@ class PointCollectionDescriptor(MeshDescriptor): units : {'degrees', 'radians'} The units of ``lats`` and ``lons`` + + history : str + The history attribute written to SCRIP files """ def __init__(self, lats, lons, collectionName, units='degrees', @@ -74,6 +75,7 @@ def __init__(self, lats, lons, collectionName, units='degrees', 'attrs': {'units': units}}} self.dims = [outDimension] self.dimSize = [len(self.lat)] + self.history = add_history() def to_scrip(self, scripFileName): """ @@ -114,7 +116,7 @@ def to_scrip(self, scripFileName): outFile.variables['grid_corner_lon'][:] = grid_corner_lon[:] # Update history attribute of netCDF file - setattr(outFile, 'history', ' '.join(sys.argv[:])) + setattr(outFile, 'history', self.history) outFile.close() diff --git a/pyremap/descriptor/projection_grid_descriptor.py b/pyremap/descriptor/projection_grid_descriptor.py index 3abaf2f..3c6e9ee 100644 --- a/pyremap/descriptor/projection_grid_descriptor.py +++ b/pyremap/descriptor/projection_grid_descriptor.py @@ -9,8 +9,6 @@ # distributed with this code, or at # https://raw.githubusercontent.com/MPAS-Dev/pyremap/main/LICENSE -import sys - import netCDF4 import numpy import pyproj @@ -18,6 +16,7 @@ from pyremap.descriptor.mesh_descriptor import MeshDescriptor from pyremap.descriptor.utility import ( + add_history, create_scrip, interp_extrap_corner, unwrap_corners, @@ -58,6 +57,9 @@ class ProjectionGridDescriptor(MeshDescriptor): yVarName : str The name of the y variable + + history : str + The history attribute written to SCRIP files """ def __init__(self, projection, meshName=None): """ @@ -84,6 +86,7 @@ def __init__(self, projection, meshName=None): self.history = None self.xVarName = None self.yVarName = None + self.history = None @classmethod def read(cls, projection, fileName, meshName=None, xVarName='x', @@ -130,12 +133,7 @@ def read(cls, projection, fileName, meshName=None, xVarName='x', descriptor.xCorner = interp_extrap_corner(descriptor.x) descriptor.yCorner = interp_extrap_corner(descriptor.y) - # Update history attribute of netCDF file - if 'history' in ds.attrs: - descriptor.history = '\n'.join([ds.attrs['history'], - ' '.join(sys.argv[:])]) - else: - descriptor.history = sys.argv[:] + descriptor.history = add_history(ds=ds) return descriptor @classmethod @@ -172,7 +170,7 @@ def create(cls, projection, x, y, meshName): # interp/extrap corners descriptor.xCorner = interp_extrap_corner(descriptor.x) descriptor.yCorner = interp_extrap_corner(descriptor.y) - descriptor.history = sys.argv[:] + descriptor.history = add_history() return descriptor def to_scrip(self, scripFileName): diff --git a/pyremap/descriptor/utility.py b/pyremap/descriptor/utility.py index 5a277d7..af481ba 100644 --- a/pyremap/descriptor/utility.py +++ b/pyremap/descriptor/utility.py @@ -9,6 +9,8 @@ # distributed with this code, or at # https://raw.githubusercontent.com/MPAS-Dev/pyremap/main/LICENSE +import sys + import numpy @@ -107,6 +109,17 @@ def unwrap_corners(inField): def round_res(res): - """Round the resoltuion to a reasonable number for grid names""" + """Round the resolution to a reasonable number for grid names""" rounded = numpy.round(res * 1000.) / 1000. return '{}'.format(rounded) + + +def add_history(ds=None): + """Get the history attribute, possibly adding it to existing history""" + history = ' '.join(sys.argv[:]) + if ds is not None and 'history' in ds.attrs: + prev_hist = ds.attrs['history'] + if isinstance(prev_hist, numpy.ndarray): + prev_hist = '\n'.join(prev_hist) + history = '\n'.join([prev_hist, history]) + return history From 608cadb6c0567beeecc8ffa2ae008592cb3095c1 Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Sat, 16 Sep 2023 22:21:07 +0200 Subject: [PATCH 6/6] Update to v1.1.0 --- ci/recipe/meta.yaml | 2 +- docs/versions.rst | 3 +++ pyremap/__init__.py | 2 +- 3 files changed, 5 insertions(+), 2 deletions(-) diff --git a/ci/recipe/meta.yaml b/ci/recipe/meta.yaml index 334ba3e..7ead132 100644 --- a/ci/recipe/meta.yaml +++ b/ci/recipe/meta.yaml @@ -1,5 +1,5 @@ {% set name = "pyremap" %} -{% set version = "1.0.1" %} +{% set version = "1.1.0" %} package: name: {{ name|lower }} diff --git a/docs/versions.rst b/docs/versions.rst index fb84800..9103ac1 100644 --- a/docs/versions.rst +++ b/docs/versions.rst @@ -18,6 +18,7 @@ Documentation On GitHub `v0.0.16`_ `0.0.16`_ `v1.0.0`_ `1.0.0`_ `v1.0.1`_ `1.0.1`_ +`v1.1.0`_ `1.1.0`_ ================ =============== .. _`stable`: ../stable/index.html @@ -34,6 +35,7 @@ Documentation On GitHub .. _`v0.0.16`: ../0.0.16/index.html .. _`v1.0.0`: ../1.0.0/index.html .. _`v1.0.1`: ../1.0.1/index.html +.. _`v1.1.0`: ../1.1.0/index.html .. _`main`: https://github.com/MPAS-Dev/pyremap/tree/main .. _`0.0.6`: https://github.com/MPAS-Dev/pyremap/tree/0.0.6 .. _`0.0.7`: https://github.com/MPAS-Dev/pyremap/tree/0.0.7 @@ -48,3 +50,4 @@ Documentation On GitHub .. _`0.0.16`: https://github.com/MPAS-Dev/pyremap/tree/0.0.16 .. _`1.0.0`: https://github.com/MPAS-Dev/pyremap/tree/1.0.0 .. _`1.0.1`: https://github.com/MPAS-Dev/pyremap/tree/1.0.1 +.. _`1.1.0`: https://github.com/MPAS-Dev/pyremap/tree/1.1.0 diff --git a/pyremap/__init__.py b/pyremap/__init__.py index e63bdd0..b7331e0 100644 --- a/pyremap/__init__.py +++ b/pyremap/__init__.py @@ -12,5 +12,5 @@ from pyremap.polar import get_polar_descriptor, get_polar_descriptor_from_file from pyremap.remapper import Remapper -__version_info__ = (1, 0, 1) +__version_info__ = (1, 1, 0) __version__ = '.'.join(str(vi) for vi in __version_info__)