Merge pull request #56 from xylar/add-mpas-vertex-mesh-descriptor

Add `MpasCellMeshDescriptor` and `MpasVertexMeshDescriptor`

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 }}

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

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,32 +32,32 @@
 # 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)
 
 # conservative remapping with 4 MPI tasks (using mpirun)
 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)

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']

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)

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)
 

pyremap/__init__.py

@@ -1,14 +1,16 @@
 from pyremap.descriptor import (
     LatLon2DGridDescriptor,
     LatLonGridDescriptor,
+    MpasCellMeshDescriptor,
     MpasEdgeMeshDescriptor,
     MpasMeshDescriptor,
+    MpasVertexMeshDescriptor,
     PointCollectionDescriptor,
     ProjectionGridDescriptor,
     get_lat_lon_descriptor,
 )
 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__)

pyremap/descriptor/__init__.py

@@ -17,8 +17,12 @@
     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.mpas_vertex_mesh_descriptor import (
+    MpasVertexMeshDescriptor,
+)
 from pyremap.descriptor.point_collection_descriptor import (
     PointCollectionDescriptor,
 )

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):

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