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

Jump to bottom

Add support for making mapping files with mbtempest #63

Open
wants to merge 7 commits into
base: main
Choose a base branch
from
Open

Add support for making mapping files with mbtempest #63

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
d98cd0d
Add moab as a dependency
xylar Jul 21, 2024
abe1761
Add moab function for making mapping files
xylar Jul 21, 2024
4797bd2
Update examples to use moab
xylar Jul 21, 2024
4eed7d0
Update API docs
xylar Jul 21, 2024
98778ff
Add 'NETCDF3_64BIT_OFFSET' and 'NETCDF3_64BIT_DATA' formats
xylar Jul 21, 2024
bffb9ea
Remove `to_esmf()` method from descriptors
xylar Jul 21, 2024
a809b99
Change default format to 'NETCDF3_64BIT_OFFSET'
xylar Jul 21, 2024
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
1 change: 1 addition & 0 deletions dev-spec.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,7 @@
python >=3.8
dask
esmf
moab=*=*_tempest_*
nco >=4.8.1
netcdf4
numpy
Expand Down
5 changes: 5 additions & 0 deletions docs/api.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -53,6 +53,11 @@ Remapping
:toctree: generated/

Remapper
Remapper.esmf_build_map
Remapper.moab_build_map
Remapper.remap_file
Remapper.remap



Polar projections
Expand Down
4 changes: 3 additions & 1 deletion examples/make_mpas_to_Antarctic_stereo_mapping.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -33,18 +33,20 @@
inGridFileName = 'ocean.QU.240km.151209.nc'

inDescriptor = MpasCellMeshDescriptor(inGridFileName, inGridName)
inDescriptor.format = 'NETCDF3_64BIT'

# 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
outDescriptor.format = 'NETCDF3_64BIT'

mappingFileName = f'map_{inGridName}_to_{outGridName}_bilinear.nc'

remapper = Remapper(inDescriptor, outDescriptor, mappingFileName)

# conservative remapping with 4 MPI tasks (using mpirun)
remapper.esmf_build_map(method='bilinear', mpi_tasks=4)
remapper.moab_build_map(method='bilinear', mpi_tasks=4)

# select the SST at the initial time as an example data set
srcFileName = f'temp_{inGridName}.nc'
Expand Down
4 changes: 3 additions & 1 deletion examples/make_mpas_to_lat_lon_mapping.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -32,17 +32,19 @@
inGridFileName = 'ocean.QU.240km.151209.nc'

inDescriptor = MpasCellMeshDescriptor(inGridFileName, inGridName)
inDescriptor.format = 'NETCDF3_64BIT'

# 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
outDescriptor.format = 'NETCDF3_64BIT'

mappingFileName = f'map_{inGridName}_to_{outGridName}_bilinear.nc'


remapper = Remapper(inDescriptor, outDescriptor, mappingFileName)

remapper.esmf_build_map(method='bilinear', mpi_tasks=1)
remapper.moab_build_map(method='bilinear', mpi_tasks=1)

outFileName = f'temp_{outGridName}.nc'
ds = xarray.open_dataset(inGridFileName)
Expand Down
6 changes: 4 additions & 2 deletions examples/make_mpas_vertex_to_Antarctic_stereo_mapping.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -34,19 +34,21 @@
inGridFileName = 'ocean.QU.240km.151209.nc'

inDescriptor = MpasVertexMeshDescriptor(inGridFileName, inGridName)
inDescriptor.format = 'NETCDF3_64BIT'

# 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
outDescriptor.format = 'NETCDF3_64BIT'

mappingFileName = f'map_{inGridName}_vertex_to_{outGridName}_conserve.nc'

remapper = Remapper(inDescriptor, outDescriptor, mappingFileName)

# conservative remapping with 4 MPI tasks (using mpirun)
remapper.esmf_build_map(method='conserve', mpi_tasks=4,
parallel_exec='mpirun', include_logs=True)
remapper.moab_build_map(method='conserve', mpi_tasks=4,
parallel_exec='mpirun')

# select the SST at the initial time as an example data set
srcFileName = f'vertex_id_{inGridName}.nc'
Expand Down
4 changes: 3 additions & 1 deletion examples/remap_stereographic.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -48,6 +48,7 @@
projection = get_antarctic_stereographic_projection()

inDescriptor = ProjectionGridDescriptor.create(projection, x, y, inMeshName)
inDescriptor.format = 'NETCDF3_64BIT'

outRes = args.resolution * 1e3

Expand All @@ -62,12 +63,13 @@

outDescriptor = ProjectionGridDescriptor.create(projection, xOut, yOut,
outMeshName)
outDescriptor.format = 'NETCDF3_64BIT'

mappingFileName = f'map_{inMeshName}_to_{outMeshName}_{args.method}.nc'

remapper = Remapper(inDescriptor, outDescriptor, mappingFileName)

remapper.esmf_build_map(method=args.method, mpi_tasks=args.mpiTasks)
remapper.moab_build_map(method=args.method, mpi_tasks=args.mpiTasks)

dsOut = remapper.remap(dsIn, renormalizationThreshold=0.01)
dsOut.to_netcdf(args.outFileName)
62 changes: 0 additions & 62 deletions pyremap/descriptor/lat_lon_grid_descriptor.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -244,68 +244,6 @@ def to_scrip(self, scripFileName, expandDist=None, expandFactor=None):

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
"""
nLat = len(self.lat)
nLon = len(self.lon)

(Lon, Lat) = numpy.meshgrid(self.lon, self.lat)
(LonCorner, LatCorner) = numpy.meshgrid(self.lonCorner, self.latCorner)
(XIndices, YIndices) = numpy.meshgrid(numpy.arange(nLon + 1),
numpy.arange(nLat + 1))

elementCount = nLon * nLat
nodeCount = (nLon + 1) * (nLat + 1)
coordDim = 2
maxNodePElement = 4

nodeCoords = numpy.zeros((nodeCount, coordDim))
nodeCoords[:, 0] = LonCorner.flat
nodeCoords[:, 1] = LatCorner.flat

centerCoords = numpy.zeros((elementCount, coordDim))
centerCoords[:, 0] = Lon.flat
centerCoords[:, 1] = Lat.flat

elementConn = numpy.zeros((elementCount, maxNodePElement), dtype=int)
elementConn[:, 0] = (XIndices[0:-1, 0:-1].ravel() +
(nLon + 1) * YIndices[0:-1, 0:-1].ravel())
elementConn[:, 1] = (XIndices[0:-1, 1:].ravel() +
(nLon + 1) * YIndices[0:-1, 1:].ravel())
elementConn[:, 2] = (XIndices[1:, 1:].ravel() +
(nLon + 1) * YIndices[1:, 1:].ravel())
elementConn[:, 3] = (XIndices[1:, 0:-1].ravel() +
(nLon + 1) * YIndices[1:, 0:-1].ravel())

ds_out = xarray.Dataset()
ds_out['nodeCoords'] = (('nodeCount', 'coordDim'), nodeCoords)
ds_out.nodeCoords.attrs['units'] = self.units
ds_out['centerCoords'] = \
(('elementCount', 'coordDim'), centerCoords)
ds_out.centerCoords.attrs['units'] = self.units
ds_out['elementConn'] = \
(('elementCount', 'maxNodePElement'), elementConn + 1)
ds_out.elementConn.attrs['start_index'] = 1
ds_out.elementConn.attrs['_FillValue'] = -1
ds_out['numElementConn'] = \
(('elementCount',), maxNodePElement * numpy.ones(elementCount,
dtype=int))

ds_out['origGridDims'] = (('origGridRank',), [nLon, nLat])

ds_out.attrs['gridType'] = 'unstructured mesh'
ds_out.attrs['version'] = '0.9'

ds_out.to_netcdf(esmfFileName, format=self.format,
engine=self.engine)

def _set_coords(self, latVarName, lonVarName, latDimName,
lonDimName):
"""
Expand Down
17 changes: 2 additions & 15 deletions pyremap/descriptor/mesh_descriptor.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -33,7 +33,7 @@ class MeshDescriptor:
coordinate in the mesh or grid

format : {'NETCDF4', 'NETCDF4_CLASSIC', 'NETCDF3_64BIT',
'NETCDF3_CLASSIC'}
'NETCDF3_64BIT_OFFSET', 'NETCDF3_64BIT_DATA', 'NETCDF3_CLASSIC'}
The NetCDF file format to use. Default is ``'NETCDF4'``

engine : {'netcdf4', 'scipy', 'h5netcdf'}
Expand All @@ -58,7 +58,7 @@ def __init__(self, meshName=None, regional=None):
self.dims = None
self.dimSizes = None
self.coords = None
self.format = 'NETCDF4'
self.format = 'NETCDF3_64BIT_OFFSET'
self.engine = None

def to_scrip(self, scripFileName, expandDist=None, expandFactor=None):
Expand All @@ -81,16 +81,3 @@ def to_scrip(self, scripFileName, expandDist=None, expandFactor=None):
"""
raise NotImplementedError(
'to_scrip is not implemented for this descriptor')

def to_esmf(self, esmfFileName):
"""
Subclasses should overload this method to write an ESMF mesh file based
on the mesh.

Parameters
----------
esmfFileName : str
The path to which the ESMF mesh file should be written
"""
raise NotImplementedError(
'to_esmf is not implemented for this descriptor')
43 changes: 0 additions & 43 deletions pyremap/descriptor/mpas_cell_mesh_descriptor.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -167,46 +167,3 @@ def to_scrip(self, scripFileName, expandDist=None, expandFactor=None):

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)
63 changes: 0 additions & 63 deletions pyremap/descriptor/mpas_edge_mesh_descriptor.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -161,66 +161,3 @@ def to_scrip(self, scripFileName, expandDist=None, expandFactor=None):

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 xr.open_dataset(self.fileName) as ds:

nCells = ds.sizes['nCells']
nodeCount = nCells + ds.sizes['nVertices']
elementCount = ds.sizes['nEdges']
coordDim = 2

nodeCoords = np.zeros((nodeCount, coordDim))
nodeCoords[0:nCells, 0] = ds.lonCell.values
nodeCoords[0:nCells, 1] = ds.latCell.values
nodeCoords[nCells:, 0] = ds.lonVertex.values
nodeCoords[nCells:, 1] = ds.latVertex.values
centerCoords = np.zeros((elementCount, coordDim))
centerCoords[:, 0] = ds.lonEdge.values
centerCoords[:, 1] = ds.latEdge.values

elementConn = -1 * np.ones((elementCount, 4), dtype=int)

verticesOnEdge = ds.verticesOnEdge.values
cellsOnEdge = ds.cellsOnEdge.values

elementConn[:, 0] = cellsOnEdge[:, 0]
elementConn[:, 1] = nCells + verticesOnEdge[:, 0]
elementConn[:, 2] = cellsOnEdge[:, 1]
elementConn[:, 3] = nCells + verticesOnEdge[:, 1]

# invalid value is -1, not 0
elementConn[elementConn == 0] = -1
for vertex in range(3):
# swap -1 value until they're at the end
mask = elementConn[:, vertex] == -1
elementConn[mask, vertex] = elementConn[mask, vertex + 1]
elementConn[mask, vertex + 1] = -1

numElementConn = np.sum((elementConn != -1), axis=1).astype(np.byte)

ds_out = xr.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',), numElementConn)
ds_out.attrs['gridType'] = 'unstructured mesh'
ds_out.attrs['version'] = '0.9'

ds_out.to_netcdf(esmfFileName, format=self.format,
engine=self.engine)
52 changes: 0 additions & 52 deletions pyremap/descriptor/point_collection_descriptor.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -11,7 +11,6 @@

import netCDF4
import numpy
import xarray

from pyremap.descriptor.mesh_descriptor import MeshDescriptor
from pyremap.descriptor.utility import add_history, create_scrip
Expand Down Expand Up @@ -127,54 +126,3 @@ def to_scrip(self, scripFileName, expandDist=None, expandFactor=None):
setattr(outFile, 'history', self.history)

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
"""
nPoints = len(self.lat)

elementCount = nPoints
nodeCount = 3 * nPoints
coordDim = 2
maxNodePElement = 3

nodeCoords = numpy.zeros((nodeCount, coordDim))
# just repeat the center lat and lon
for iVertex in range(maxNodePElement):
nodeCoords[iVertex::maxNodePElement, 0] = self.lon
nodeCoords[iVertex::maxNodePElement, 1] = self.lat

centerCoords = numpy.zeros((elementCount, coordDim))
centerCoords[:, 0] = self.lon
centerCoords[:, 1] = self.lat

elementConn = numpy.zeros((elementCount, maxNodePElement), dtype=int)
elementConn[:, 0] = maxNodePElement * numpy.arange(nPoints)
elementConn[:, 1] = maxNodePElement * numpy.arange(nPoints) + 1
elementConn[:, 2] = maxNodePElement * numpy.arange(nPoints) + 2

ds_out = xarray.Dataset()
ds_out['nodeCoords'] = (('nodeCount', 'coordDim'), nodeCoords)
ds_out.nodeCoords.attrs['units'] = self.units
ds_out['centerCoords'] = \
(('elementCount', 'coordDim'), centerCoords)
ds_out.centerCoords.attrs['units'] = self.units
ds_out['elementConn'] = \
(('elementCount', 'maxNodePElement'), elementConn + 1)
ds_out.elementConn.attrs['start_index'] = 1
ds_out.elementConn.attrs['_FillValue'] = -1
ds_out['numElementConn'] = \
(('elementCount',), maxNodePElement * numpy.ones(elementCount,
dtype=int))

ds_out.attrs['gridType'] = 'unstructured mesh'
ds_out.attrs['version'] = '0.9'

ds_out.to_netcdf(esmfFileName, format=self.format,
engine=self.engine)
Loading
Loading