From d98cd0d71e52bbcc30f77c2ed125c8261a1c5cd6 Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Sun, 21 Jul 2024 14:28:08 +0200 Subject: [PATCH 1/7] Add moab as a dependency --- dev-spec.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/dev-spec.txt b/dev-spec.txt index f7a2deb..058a57f 100644 --- a/dev-spec.txt +++ b/dev-spec.txt @@ -5,6 +5,7 @@ python >=3.8 dask esmf +moab=*=*_tempest_* nco >=4.8.1 netcdf4 numpy From abe1761bf9e1717b7f7009cefda07a0652bdba28 Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Sun, 21 Jul 2024 14:31:41 +0200 Subject: [PATCH 2/7] Add moab function for making mapping files --- pyremap/remapper.py | 235 +++++++++++++++++++++++++++++++++++++------- 1 file changed, 202 insertions(+), 33 deletions(-) diff --git a/pyremap/remapper.py b/pyremap/remapper.py index 813b4f5..d8e8215 100644 --- a/pyremap/remapper.py +++ b/pyremap/remapper.py @@ -156,6 +156,63 @@ def esmf_build_map(self, method='bilinear', logger=None, mpi_tasks=1, expand_dist=expand_dist, expand_factor=expand_factor) + def moab_build_map(self, method='bilinear', logger=None, mpi_tasks=1, + tempdir=None, parallel_exec=None, expand_dist=None, + expand_factor=None): + """ + Use ``mbtempest`` to construct a mapping file used for interpolation + between the source and destination grids. + + Parameters + ---------- + method : {'bilinear', 'conserve'}, optional + The method of interpolation used, see documentation for + ``mbtempest`` for details. + + logger : ``logging.Logger``, optional + A logger to which ncclimo output should be redirected + + mpi_tasks : int, optional + The number of MPI tasks (a number > 1 implies that ``mbtempest`` + will be called with the given parallel executable) + + tempdir : str, optional + A temporary directory. By default, a temporary directory is + created, typically in ``/tmp`` but on some systems such as compute + nodes this may not be visible to all processors in the subsequent + ``ESMF_RegridWeightGen`` call + + parallel_exec : {'srun', 'mpirun'}, optional + The name of the parallel executable to use to launch ``mbtempest``. + By default, 'mpirun' from the conda environment is used + + expand_dist : float or numpy.ndarray, optional + A distance in meters to expand each cell of the destination mesh + outward from the center. This can be used to smooth fields on + the destination grid. If a ``numpy.ndarray``, one value per cell + on the destination mesh. + + expand_factor : float or numpy.ndarray, optional + A factor by which to expand each cell of the destination mesh + outward from the center. This can be used to smooth fields on + the destination grid. If a ``numpy.ndarray``, one value per cell + on the destination mesh. + + Raises + ------ + OSError + If ``mbtempest`` is not in the system path. + + ValueError + If the type of mesh used in ``sourceDescriptor`` and/or + ``destinationDescriptor`` are not compatible with other arguments + """ + self._moab_build_map(method=method, logger=logger, + mpi_tasks=mpi_tasks, + tempdir=tempdir, parallel_exec=parallel_exec, + expand_dist=expand_dist, + expand_factor=expand_factor) + def build_mapping_file(self, method='bilinear', additionalArgs=None, logger=None, mpiTasks=1, tempdir=None, esmf_path=None, @@ -230,7 +287,7 @@ def build_mapping_file(self, method='bilinear', # Xylar Asay-Davis warn('build_mapping_file() is deprecated. Use esmf_build_map() or ' - 'mbtempest_build_map() instead', DeprecationWarning, stacklevel=2) + 'moab_build_map() instead', DeprecationWarning, stacklevel=2) self._esmf_build_map(method=method, logger=logger, mpi_tasks=mpiTasks, tempdir=tempdir, parallel_exec=esmf_parallel_exec, @@ -759,42 +816,13 @@ def _esmf_build_map(self, method, logger, mpi_tasks, tempdir, if additional_args is not None: args.extend(additional_args) - if logger is None: - _print_running(args, fn=print) - # make sure any output is flushed before we add output from the - # subprocess - sys.stdout.flush() - sys.stderr.flush() - - # throw out the standard output from ESMF_RegridWeightGen, as it's - # rather verbose but keep stderr - with open(os.devnull, 'wb') as DEVNULL: - subprocess.check_call(args, stdout=DEVNULL) - - else: - _print_running(args, fn=logger.info) - for handler in logger.handlers: - handler.flush() - - process = subprocess.Popen(args, stdout=subprocess.PIPE, - stderr=subprocess.PIPE) - _, stderr = process.communicate() - - # throw out the standard output from ESMF_RegridWeightGen, as it's - # rather verbose but keep stderr - if stderr: - stderr = stderr.decode('utf-8') - for line in stderr.split('\n'): - logger.error(line) - - if process.returncode != 0: - raise subprocess.CalledProcessError(process.returncode, - ' '.join(args)) + _run_subprocess(args, logger) if tempobj is not None: tempobj.cleanup() - def _validate_descriptors(self, method, expand_factor, expand_dist): + def _validate_descriptors(self, method, expand_factor, expand_dist, + netcdf3=False): if isinstance(self.destinationDescriptor, PointCollectionDescriptor) and \ @@ -807,6 +835,17 @@ def _validate_descriptors(self, method, expand_factor, expand_dist): raise ValueError('expandFactor and expandDist only have an impact ' 'with method=conserve') + if netcdf3: + allowed_formats = ['NETCDF3_64BIT', 'NETCDF3_CLASSIC'] + if self.sourceDescriptor.format not in allowed_formats: + raise ValueError(f'sourceDescriptor must be in one of the ' + f'formats {allowed_formats} but it is ' + f'{self.sourceDescriptor.format}') + if self.destinationDescriptor.format not in allowed_formats: + raise ValueError(f'destinationDescriptor must be in one of ' + f'the formats {allowed_formats} but it is ' + f'{self.destinationDescriptor.format}') + def _get_rwg_path(self, esmf_path): if esmf_path is not None: @@ -872,6 +911,103 @@ def _get_esmf_regional_args(self): return regional_args + def _moab_build_map(self, method, logger, mpi_tasks, tempdir, + parallel_exec, expand_dist, expand_factor): + """ + Use ``mbtempest`` to construct a mapping file used for interpolation + between the source and destination grids. + """ + + if self.mappingFileName is None or \ + os.path.exists(self.mappingFileName): + # a valid weight file already exists, so nothing to do + return + + self._validate_descriptors(method, expand_factor, expand_dist, + netcdf3=True) + + mbtempest_path = self._get_mbtempest_path() + + # Write source and destination SCRIP files in temporary locations + if tempdir is None: + tempobj = TemporaryDirectory() + tempdir = tempobj.name + else: + tempobj = None + + src_filename = f'{tempdir}/src_mesh.nc' + dst_filename = f'{tempdir}/dst_mesh.nc' + src_descriptor = self.sourceDescriptor + dst_descriptor = self.destinationDescriptor + + src_descriptor.to_scrip(src_filename) + + dst_descriptor.to_scrip(dst_filename, + expandDist=expand_dist, + expandFactor=expand_factor) + + parallel_args = self._get_parallel_args(parallel_exec, mpi_tasks, + conda_package='moab') + regional_args = self._get_mbtempest_regional_args() + + intx_filename = \ + f'moab_intx_{src_descriptor.meshName}_to_' \ + f'{dst_descriptor.meshName}.h5m' + + intx_args = parallel_args + [ + mbtempest_path, + '--type', '5', + '--load', src_filename, + '--load', dst_filename, + '--intx', intx_filename + ] + regional_args + + _run_subprocess(intx_args, logger) + + fvmethod = { + 'conserve': 'none', + 'bilinear': 'bilin'} + + map_args = parallel_args + [ + mbtempest_path, + '--type', '5', + '--load', src_filename, + '--load', dst_filename, + '--intx', intx_filename, + '--weights', + '--method', 'fv', + '--method', 'fv', + '--file', self.mappingFileName, + '--order', '1', + '--order', '1', + '--fvmethod', fvmethod[method] + ] + regional_args + + _run_subprocess(map_args, logger) + + if tempobj is not None: + tempobj.cleanup() + + def _get_mbtempest_path(self): + + mbtempest_path = find_executable('mbtempest') + + if mbtempest_path is None: + raise OSError('mbtempest not found. Make sure moab package is ' + 'installed with tempest support: \n' + 'conda install "moab=*=*_tempest_*"\n') + return mbtempest_path + + def _get_mbtempest_regional_args(self): + + regional_args = [] + + if self.sourceDescriptor.regional or \ + self.destinationDescriptor.regional: + regional_args.append('--rrmgrids') + + return regional_args + def _print_running(args, fn): print_args = [] @@ -880,3 +1016,36 @@ def _print_running(args, fn): arg = f'"{arg}"' print_args.append(arg) fn(f'running: {" ".join(print_args)}') + + +def _run_subprocess(args, logger): + if logger is None: + log_err_fn = print + log_out_fn = print + _print_running(args, fn=log_out_fn) + sys.stdout.flush() + sys.stderr.flush() + else: + log_err_fn = logger.error + log_out_fn = logger.info + _print_running(args, fn=log_out_fn) + for handler in logger.handlers: + handler.flush() + + with subprocess.Popen(args, stdout=subprocess.PIPE, + stderr=subprocess.PIPE) as process: + stdout, stderr = process.communicate() + + if stdout: + stdout = stdout.decode('utf-8') + for line in stdout.split('\n'): + log_out_fn(line) + + if stderr: + stderr = stderr.decode('utf-8') + for line in stderr.split('\n'): + log_err_fn(line) + + if process.returncode != 0: + raise subprocess.CalledProcessError(process.returncode, + ' '.join(args)) From 4797bd23c4699117277b04bc68c651134cc5824b Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Sun, 21 Jul 2024 14:32:02 +0200 Subject: [PATCH 3/7] Update examples to use moab --- examples/make_mpas_to_Antarctic_stereo_mapping.py | 4 +++- examples/make_mpas_to_lat_lon_mapping.py | 4 +++- examples/make_mpas_vertex_to_Antarctic_stereo_mapping.py | 6 ++++-- examples/remap_stereographic.py | 4 +++- 4 files changed, 13 insertions(+), 5 deletions(-) diff --git a/examples/make_mpas_to_Antarctic_stereo_mapping.py b/examples/make_mpas_to_Antarctic_stereo_mapping.py index d9262c1..1440e0c 100755 --- a/examples/make_mpas_to_Antarctic_stereo_mapping.py +++ b/examples/make_mpas_to_Antarctic_stereo_mapping.py @@ -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' diff --git a/examples/make_mpas_to_lat_lon_mapping.py b/examples/make_mpas_to_lat_lon_mapping.py index adf5f49..dec9405 100755 --- a/examples/make_mpas_to_lat_lon_mapping.py +++ b/examples/make_mpas_to_lat_lon_mapping.py @@ -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) diff --git a/examples/make_mpas_vertex_to_Antarctic_stereo_mapping.py b/examples/make_mpas_vertex_to_Antarctic_stereo_mapping.py index 9c34952..7e0173e 100755 --- a/examples/make_mpas_vertex_to_Antarctic_stereo_mapping.py +++ b/examples/make_mpas_vertex_to_Antarctic_stereo_mapping.py @@ -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' diff --git a/examples/remap_stereographic.py b/examples/remap_stereographic.py index 5bcd5dc..7c02fe7 100755 --- a/examples/remap_stereographic.py +++ b/examples/remap_stereographic.py @@ -48,6 +48,7 @@ projection = get_antarctic_stereographic_projection() inDescriptor = ProjectionGridDescriptor.create(projection, x, y, inMeshName) +inDescriptor.format = 'NETCDF3_64BIT' outRes = args.resolution * 1e3 @@ -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) From 4eed7d0ee5722d53e70f119280fb1032aca5338f Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Sun, 21 Jul 2024 14:35:30 +0200 Subject: [PATCH 4/7] Update API docs --- docs/api.rst | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/docs/api.rst b/docs/api.rst index 7a86c10..9ddca03 100644 --- a/docs/api.rst +++ b/docs/api.rst @@ -53,6 +53,11 @@ Remapping :toctree: generated/ Remapper + Remapper.esmf_build_map + Remapper.moab_build_map + Remapper.remap_file + Remapper.remap + Polar projections From 98778ff694a74660b1b94aff379203a770201223 Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Sun, 21 Jul 2024 17:37:23 +0200 Subject: [PATCH 5/7] Add 'NETCDF3_64BIT_OFFSET' and 'NETCDF3_64BIT_DATA' formats The latter is needed for large input grids to MOAB. --- pyremap/descriptor/mesh_descriptor.py | 2 +- pyremap/remapper.py | 3 ++- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/pyremap/descriptor/mesh_descriptor.py b/pyremap/descriptor/mesh_descriptor.py index d24d696..e54b019 100644 --- a/pyremap/descriptor/mesh_descriptor.py +++ b/pyremap/descriptor/mesh_descriptor.py @@ -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'} diff --git a/pyremap/remapper.py b/pyremap/remapper.py index d8e8215..4defb13 100644 --- a/pyremap/remapper.py +++ b/pyremap/remapper.py @@ -836,7 +836,8 @@ def _validate_descriptors(self, method, expand_factor, expand_dist, 'with method=conserve') if netcdf3: - allowed_formats = ['NETCDF3_64BIT', 'NETCDF3_CLASSIC'] + allowed_formats = ['NETCDF3_64BIT', 'NETCDF3_64BIT_OFFSET', + 'NETCDF3_64BIT_DATA', 'NETCDF3_CLASSIC'] if self.sourceDescriptor.format not in allowed_formats: raise ValueError(f'sourceDescriptor must be in one of the ' f'formats {allowed_formats} but it is ' From bffb9eaae128fb1290584e895a09d8fb8a449cd4 Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Sun, 21 Jul 2024 17:42:37 +0200 Subject: [PATCH 6/7] Remove `to_esmf()` method from descriptors MOAB cannot read these grids so we don't want to put extra effort into supporting them. --- pyremap/descriptor/lat_lon_grid_descriptor.py | 62 ------------------ pyremap/descriptor/mesh_descriptor.py | 13 ---- .../descriptor/mpas_cell_mesh_descriptor.py | 43 ------------ .../descriptor/mpas_edge_mesh_descriptor.py | 63 ------------------ .../descriptor/point_collection_descriptor.py | 52 --------------- .../descriptor/projection_grid_descriptor.py | 65 ------------------- 6 files changed, 298 deletions(-) diff --git a/pyremap/descriptor/lat_lon_grid_descriptor.py b/pyremap/descriptor/lat_lon_grid_descriptor.py index 8b55744..b5ce15a 100644 --- a/pyremap/descriptor/lat_lon_grid_descriptor.py +++ b/pyremap/descriptor/lat_lon_grid_descriptor.py @@ -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): """ diff --git a/pyremap/descriptor/mesh_descriptor.py b/pyremap/descriptor/mesh_descriptor.py index e54b019..d145c5a 100644 --- a/pyremap/descriptor/mesh_descriptor.py +++ b/pyremap/descriptor/mesh_descriptor.py @@ -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') diff --git a/pyremap/descriptor/mpas_cell_mesh_descriptor.py b/pyremap/descriptor/mpas_cell_mesh_descriptor.py index 7c88de3..4804d62 100644 --- a/pyremap/descriptor/mpas_cell_mesh_descriptor.py +++ b/pyremap/descriptor/mpas_cell_mesh_descriptor.py @@ -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) diff --git a/pyremap/descriptor/mpas_edge_mesh_descriptor.py b/pyremap/descriptor/mpas_edge_mesh_descriptor.py index 75887dd..72d0d00 100644 --- a/pyremap/descriptor/mpas_edge_mesh_descriptor.py +++ b/pyremap/descriptor/mpas_edge_mesh_descriptor.py @@ -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) diff --git a/pyremap/descriptor/point_collection_descriptor.py b/pyremap/descriptor/point_collection_descriptor.py index 730de3a..2107601 100644 --- a/pyremap/descriptor/point_collection_descriptor.py +++ b/pyremap/descriptor/point_collection_descriptor.py @@ -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 @@ -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) diff --git a/pyremap/descriptor/projection_grid_descriptor.py b/pyremap/descriptor/projection_grid_descriptor.py index a58fee6..7b78b81 100644 --- a/pyremap/descriptor/projection_grid_descriptor.py +++ b/pyremap/descriptor/projection_grid_descriptor.py @@ -222,71 +222,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 - """ - nx = len(self.x) - ny = len(self.y) - - (X, Y) = numpy.meshgrid(self.x, self.y) - (XCorner, YCorner) = numpy.meshgrid(self.xCorner, self.yCorner) - (XIndices, YIndices) = numpy.meshgrid(numpy.arange(nx + 1), - numpy.arange(ny + 1)) - - (Lat, Lon) = self.project_to_lat_lon(X, Y) - (LatCorner, LonCorner) = self.project_to_lat_lon(XCorner, YCorner) - - elementCount = nx * ny - nodeCount = (nx + 1) * (ny + 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() + - (nx + 1) * YIndices[0:-1, 0:-1].ravel()) - elementConn[:, 1] = (XIndices[0:-1, 1:].ravel() + - (nx + 1) * YIndices[0:-1, 1:].ravel()) - elementConn[:, 2] = (XIndices[1:, 1:].ravel() + - (nx + 1) * YIndices[1:, 1:].ravel()) - elementConn[:, 3] = (XIndices[1:, 0:-1].ravel() + - (nx + 1) * YIndices[1:, 0:-1].ravel()) - - ds_out = xarray.Dataset() - ds_out['nodeCoords'] = (('nodeCount', 'coordDim'), nodeCoords) - ds_out.nodeCoords.attrs['units'] = 'degrees' - ds_out['centerCoords'] = \ - (('elementCount', 'coordDim'), centerCoords) - ds_out.centerCoords.attrs['units'] = 'degrees' - 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',), [ny, nx]) - - ds_out.attrs['gridType'] = 'unstructured mesh' - ds_out.attrs['version'] = '0.9' - - ds_out.to_netcdf(esmfFileName, format=self.format, - engine=self.engine) - def project_to_lat_lon(self, X, Y): """ Given X and Y locations of points in a projection, returns the From a809b994f288551eda97bfecd7b594aea3714f31 Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Sun, 21 Jul 2024 17:48:37 +0200 Subject: [PATCH 7/7] Change default format to 'NETCDF3_64BIT_OFFSET' This will beetter accommodate MOAB --- pyremap/descriptor/mesh_descriptor.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyremap/descriptor/mesh_descriptor.py b/pyremap/descriptor/mesh_descriptor.py index d145c5a..1bb9065 100644 --- a/pyremap/descriptor/mesh_descriptor.py +++ b/pyremap/descriptor/mesh_descriptor.py @@ -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):