diff --git a/ci/recipe/meta.yaml b/ci/recipe/meta.yaml index 7ead132..cdd8b15 100644 --- a/ci/recipe/meta.yaml +++ b/ci/recipe/meta.yaml @@ -1,5 +1,5 @@ {% set name = "pyremap" %} -{% set version = "1.1.0" %} +{% set version = "1.2.0" %} package: name: {{ name|lower }} diff --git a/docs/versions.rst b/docs/versions.rst index 9103ac1..c6874e4 100644 --- a/docs/versions.rst +++ b/docs/versions.rst @@ -19,6 +19,7 @@ Documentation On GitHub `v1.0.0`_ `1.0.0`_ `v1.0.1`_ `1.0.1`_ `v1.1.0`_ `1.1.0`_ +`v1.2.0`_ `1.2.0`_ ================ =============== .. _`stable`: ../stable/index.html @@ -36,6 +37,7 @@ Documentation On GitHub .. _`v1.0.0`: ../1.0.0/index.html .. _`v1.0.1`: ../1.0.1/index.html .. _`v1.1.0`: ../1.1.0/index.html +.. _`v1.2.0`: ../1.2.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 @@ -51,3 +53,4 @@ Documentation On GitHub .. _`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 +.. _`1.2.0`: https://github.com/MPAS-Dev/pyremap/tree/1.2.0 diff --git a/pyremap/__init__.py b/pyremap/__init__.py index b7331e0..a488ca5 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, 1, 0) +__version_info__ = (1, 2, 0) __version__ = '.'.join(str(vi) for vi in __version_info__) diff --git a/pyremap/descriptor/lat_lon_2d_grid_descriptor.py b/pyremap/descriptor/lat_lon_2d_grid_descriptor.py index 8e8727a..4d2aaea 100644 --- a/pyremap/descriptor/lat_lon_2d_grid_descriptor.py +++ b/pyremap/descriptor/lat_lon_2d_grid_descriptor.py @@ -17,6 +17,7 @@ from pyremap.descriptor.utility import ( add_history, create_scrip, + expand_scrip, interp_extrap_corners_2d, round_res, unwrap_corners, @@ -122,7 +123,7 @@ def read(cls, fileName=None, ds=None, latVarName='lat', descriptor.history = add_history(ds=ds) return descriptor - def to_scrip(self, scripFileName): + def to_scrip(self, scripFileName, expandDist=None, expandFactor=None): """ Create a SCRIP file based on the grid. @@ -130,6 +131,13 @@ def to_scrip(self, scripFileName): ---------- scripFileName : str The path to which the SCRIP file should be written + + expandDist : float, optional + A distance in meters to expand each grid cell outward from the + center + + expandFactor : float, optional + A factor by which to expand each grid cell outward from the center """ outFile = netCDF4.Dataset(scripFileName, 'w', format=self.format) @@ -150,6 +158,9 @@ def to_scrip(self, scripFileName): outFile.variables['grid_corner_lon'][:] = \ unwrap_corners(self.lonCorner) + if expandDist is not None or expandFactor is not None: + expand_scrip(outFile, expandDist, expandFactor) + setattr(outFile, 'history', self.history) outFile.close() diff --git a/pyremap/descriptor/lat_lon_grid_descriptor.py b/pyremap/descriptor/lat_lon_grid_descriptor.py index 6fb4d2b..1511442 100644 --- a/pyremap/descriptor/lat_lon_grid_descriptor.py +++ b/pyremap/descriptor/lat_lon_grid_descriptor.py @@ -17,6 +17,7 @@ from pyremap.descriptor.utility import ( add_history, create_scrip, + expand_scrip, interp_extrap_corner, round_res, unwrap_corners, @@ -198,7 +199,7 @@ def create(cls, latCorner, lonCorner, units='degrees', meshName=None, descriptor._set_coords('lat', 'lon', 'lat', 'lon') return descriptor - def to_scrip(self, scripFileName): + def to_scrip(self, scripFileName, expandDist=None, expandFactor=None): """ Given a lat-lon grid file, create a SCRIP file based on the grid. @@ -206,6 +207,13 @@ def to_scrip(self, scripFileName): ---------- scripFileName : str The path to which the SCRIP file should be written + + expandDist : float, optional + A distance in meters to expand each grid cell outward from the + center + + expandFactor : float, optional + A factor by which to expand each grid cell outward from the center """ outFile = netCDF4.Dataset(scripFileName, 'w', format=self.format) @@ -228,6 +236,9 @@ def to_scrip(self, scripFileName): outFile.variables['grid_corner_lat'][:] = unwrap_corners(LatCorner) outFile.variables['grid_corner_lon'][:] = unwrap_corners(LonCorner) + if expandDist is not None or expandFactor is not None: + expand_scrip(outFile, expandDist, expandFactor) + setattr(outFile, 'history', self.history) outFile.close() diff --git a/pyremap/descriptor/mesh_descriptor.py b/pyremap/descriptor/mesh_descriptor.py index 7c8cbb4..fcb8349 100644 --- a/pyremap/descriptor/mesh_descriptor.py +++ b/pyremap/descriptor/mesh_descriptor.py @@ -61,7 +61,7 @@ def __init__(self, meshName=None, regional=None): self.format = 'NETCDF4' self.engine = None - def to_scrip(self, scripFileName): + def to_scrip(self, scripFileName, expandDist=None, expandFactor=None): """ Subclasses should overload this method to write a SCRIP file based on the mesh. @@ -70,6 +70,13 @@ def to_scrip(self, scripFileName): ---------- scripFileName : str The path to which the SCRIP file should be written + + expandDist : float, optional + A distance in meters to expand each grid cell outward from the + center + + expandFactor : float, optional + A factor by which to expand each grid cell outward from the center """ raise NotImplementedError( 'to_scrip 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 904f76d..185ab27 100644 --- a/pyremap/descriptor/mpas_cell_mesh_descriptor.py +++ b/pyremap/descriptor/mpas_cell_mesh_descriptor.py @@ -16,7 +16,7 @@ import xarray from pyremap.descriptor.mesh_descriptor import MeshDescriptor -from pyremap.descriptor.utility import add_history, create_scrip +from pyremap.descriptor.utility import add_history, create_scrip, expand_scrip class MpasCellMeshDescriptor(MeshDescriptor): @@ -98,7 +98,7 @@ def __init__(self, fileName, meshName=None, vertices=False): self.history = add_history(ds=ds) - def to_scrip(self, scripFileName): + def to_scrip(self, scripFileName, expandDist=None, expandFactor=None): """ Given an MPAS mesh file, create a SCRIP file based on the mesh. @@ -106,6 +106,13 @@ def to_scrip(self, scripFileName): ---------- scripFileName : str The path to which the SCRIP file should be written + + expandDist : float, optional + A distance in meters to expand each grid cell outward from the + center + + expandFactor : float, optional + A factor by which to expand each grid cell outward from the center """ if self.vertices: raise ValueError('A SCRIP file won\'t work for remapping vertices') @@ -152,6 +159,9 @@ def to_scrip(self, scripFileName): outFile.variables['grid_corner_lat'][:] = grid_corner_lat[:] outFile.variables['grid_corner_lon'][:] = grid_corner_lon[:] + if expandDist is not None or expandFactor is not None: + expand_scrip(outFile, expandDist, expandFactor) + setattr(outFile, 'history', self.history) inFile.close() diff --git a/pyremap/descriptor/mpas_edge_mesh_descriptor.py b/pyremap/descriptor/mpas_edge_mesh_descriptor.py index 0d385c6..85571ff 100644 --- a/pyremap/descriptor/mpas_edge_mesh_descriptor.py +++ b/pyremap/descriptor/mpas_edge_mesh_descriptor.py @@ -14,7 +14,7 @@ import xarray as xr from pyremap.descriptor.mesh_descriptor import MeshDescriptor -from pyremap.descriptor.utility import add_history, create_scrip +from pyremap.descriptor.utility import add_history, create_scrip, expand_scrip class MpasEdgeMeshDescriptor(MeshDescriptor): @@ -70,7 +70,7 @@ def __init__(self, fileName, meshName=None): self.history = add_history(ds=ds) - def to_scrip(self, scripFileName): + def to_scrip(self, scripFileName, expandDist=None, expandFactor=None): """ Given an MPAS mesh file, create a SCRIP file based on the mesh. @@ -78,6 +78,13 @@ def to_scrip(self, scripFileName): ---------- scripFileName : str The path to which the SCRIP file should be written + + expandDist : float, optional + A distance in meters to expand each grid cell outward from the + center + + expandFactor : float, optional + A factor by which to expand each grid cell outward from the center """ inFile = netCDF4.Dataset(self.fileName, 'r') @@ -146,6 +153,9 @@ def to_scrip(self, scripFileName): outFile.variables['grid_corner_lat'][:] = grid_corner_lat[:] outFile.variables['grid_corner_lon'][:] = grid_corner_lon[:] + if expandDist is not None or expandFactor is not None: + expand_scrip(outFile, expandDist, expandFactor) + setattr(outFile, 'history', self.history) inFile.close() diff --git a/pyremap/descriptor/mpas_vertex_mesh_descriptor.py b/pyremap/descriptor/mpas_vertex_mesh_descriptor.py index 0e6682f..d7afe91 100644 --- a/pyremap/descriptor/mpas_vertex_mesh_descriptor.py +++ b/pyremap/descriptor/mpas_vertex_mesh_descriptor.py @@ -14,7 +14,7 @@ import xarray as xr from pyremap.descriptor.mesh_descriptor import MeshDescriptor -from pyremap.descriptor.utility import add_history, create_scrip +from pyremap.descriptor.utility import add_history, create_scrip, expand_scrip class MpasVertexMeshDescriptor(MeshDescriptor): @@ -70,7 +70,7 @@ def __init__(self, fileName, meshName=None): self.history = add_history(ds=ds) - def to_scrip(self, scripFileName): + def to_scrip(self, scripFileName, expandDist=None, expandFactor=None): """ Given an MPAS mesh file, create a SCRIP file based on the mesh. @@ -78,6 +78,13 @@ def to_scrip(self, scripFileName): ---------- scripFileName : str The path to which the SCRIP file should be written + + expandDist : float, optional + A distance in meters to expand each grid cell outward from the + center + + expandFactor : float, optional + A factor by which to expand each grid cell outward from the center """ inFile = netCDF4.Dataset(self.fileName, 'r') @@ -150,6 +157,9 @@ def to_scrip(self, scripFileName): outFile.variables['grid_corner_lat'][:] = grid_corner_lat[:] outFile.variables['grid_corner_lon'][:] = grid_corner_lon[:] + if expandDist is not None or expandFactor is not None: + expand_scrip(outFile, expandDist, expandFactor) + setattr(outFile, 'history', self.history) inFile.close() diff --git a/pyremap/descriptor/point_collection_descriptor.py b/pyremap/descriptor/point_collection_descriptor.py index 679e7e3..03dc372 100644 --- a/pyremap/descriptor/point_collection_descriptor.py +++ b/pyremap/descriptor/point_collection_descriptor.py @@ -77,7 +77,7 @@ def __init__(self, lats, lons, collectionName, units='degrees', self.dimSize = [len(self.lat)] self.history = add_history() - def to_scrip(self, scripFileName): + def to_scrip(self, scripFileName, expandDist=None, expandFactor=None): """ Given an MPAS mesh file, create a SCRIP file based on the mesh. @@ -85,6 +85,13 @@ def to_scrip(self, scripFileName): ---------- scripFileName : str The path to which the SCRIP file should be written + + expandDist : float, optional + A distance in meters to expand each grid cell outward from the + center + + expandFactor : float, optional + A factor by which to expand each grid cell outward from the center """ outFile = netCDF4.Dataset(scripFileName, 'w', format=self.format) diff --git a/pyremap/descriptor/projection_grid_descriptor.py b/pyremap/descriptor/projection_grid_descriptor.py index 3c6e9ee..aef0235 100644 --- a/pyremap/descriptor/projection_grid_descriptor.py +++ b/pyremap/descriptor/projection_grid_descriptor.py @@ -18,6 +18,7 @@ from pyremap.descriptor.utility import ( add_history, create_scrip, + expand_scrip, interp_extrap_corner, unwrap_corners, ) @@ -173,7 +174,7 @@ def create(cls, projection, x, y, meshName): descriptor.history = add_history() return descriptor - def to_scrip(self, scripFileName): + def to_scrip(self, scripFileName, expandDist=None, expandFactor=None): """ Create a SCRIP file based on the grid and projection. @@ -181,6 +182,13 @@ def to_scrip(self, scripFileName): ---------- scripFileName : str The path to which the SCRIP file should be written + + expandDist : float, optional + A distance in meters to expand each grid cell outward from the + center + + expandFactor : float, optional + A factor by which to expand each grid cell outward from the center """ outFile = netCDF4.Dataset(scripFileName, 'w', format=self.format) @@ -206,6 +214,9 @@ def to_scrip(self, scripFileName): outFile.variables['grid_corner_lat'][:] = unwrap_corners(LatCorner) outFile.variables['grid_corner_lon'][:] = unwrap_corners(LonCorner) + if expandDist is not None or expandFactor is not None: + expand_scrip(outFile, expandDist, expandFactor) + setattr(outFile, 'history', self.history) outFile.close() diff --git a/pyremap/descriptor/utility.py b/pyremap/descriptor/utility.py index af481ba..0c55679 100644 --- a/pyremap/descriptor/utility.py +++ b/pyremap/descriptor/utility.py @@ -12,6 +12,8 @@ import sys import numpy +import pyproj.enums +from pyproj import Transformer def interp_extrap_corner(inField): @@ -95,6 +97,82 @@ def create_scrip(outFile, grid_size, grid_corners, grid_rank, units, outFile.meshName = meshName +def expand_scrip(outFile, expandDist, expandFactor): + """ + Given a SCRIP files, creates common variables and writes common values used + in various types of SCRIP files. + + Parameters + ---------- + outFile : file pointer + A SCRIP file opened in write mode + + expandDist : float + A distance in meters to expand each grid cell outward from the + center + + expandFactor : float + A factor by which to expand each grid cell outward from the center + """ + grid_center_lat = outFile.variables['grid_center_lat'][:] + grid_center_lon = outFile.variables['grid_center_lon'][:] + grid_corner_lat = outFile.variables['grid_corner_lat'][:] + grid_corner_lon = outFile.variables['grid_corner_lon'][:] + grid_size = len(outFile.dimensions['grid_size']) + grid_corners = len(outFile.dimensions['grid_corners']) + + radians = 'rad' in outFile.variables['grid_center_lat'].units + + print(radians) + + for index in range(grid_corners): + print(index) + print(grid_corner_lon[0, index], grid_corner_lat[0, index]) + + trans_lon_lat_to_xyz = Transformer.from_crs(4979, 4978, always_xy=True) + x_center, y_center, z_center = trans_lon_lat_to_xyz.transform( + grid_center_lon, grid_center_lat, numpy.zeros(grid_size), + radians=radians) + + x_corner, y_corner, z_corner = trans_lon_lat_to_xyz.transform( + grid_corner_lon, grid_corner_lat, + numpy.zeros((grid_size, grid_corners)), radians=radians) + + if expandFactor is None: + expandFactor = 1. + + if expandDist is None: + expandDist = 0. + + for index in range(grid_corners): + print(index) + print(x_corner[0, index], y_corner[0, index], z_corner[0, index]) + print(x_center[0], y_center[0], z_center[0]) + dx = x_corner[:, index] - x_center + dy = y_corner[:, index] - y_center + dz = z_corner[:, index] - z_center + print(dx[0], dy[0], dz[0]) + dist = numpy.sqrt(dx**2 + dy**2 + dz**2) + print(dist[0]) + factor = (expandFactor * dist + expandDist) / dist + print(factor[0]) + x_corner[:, index] = factor * dx + x_center + y_corner[:, index] = factor * dy + y_center + z_corner[:, index] = factor * dz + z_center + print(x_corner[0, index], y_corner[0, index], z_corner[0, index]) + + grid_corner_lon, grid_corner_lat, _ = trans_lon_lat_to_xyz.transform( + x_corner, y_corner, z_corner, radians=radians, + direction=pyproj.enums.TransformDirection.INVERSE) + + for index in range(grid_corners): + print(index) + print(grid_corner_lon[0, index], grid_corner_lat[0, index]) + + outFile.variables['grid_corner_lat'][:] = grid_corner_lat[:] + outFile.variables['grid_corner_lon'][:] = grid_corner_lon[:] + + def unwrap_corners(inField): """Turn a 2D array of corners into an array of rectangular mesh elements""" outField = numpy.zeros(((inField.shape[0] - 1) * (inField.shape[1] - 1), diff --git a/pyremap/remapper.py b/pyremap/remapper.py index ac02d53..f730f4e 100644 --- a/pyremap/remapper.py +++ b/pyremap/remapper.py @@ -98,7 +98,8 @@ 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, - include_logs=False): + include_logs=False, expandDist=None, + expandFactor=None): """ 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 @@ -142,6 +143,16 @@ def build_mapping_file(self, method='bilinear', # noqa: C901 include_logs : bool, optional Whether to include log files from ``ESMF_RegridWeightGen`` + expandDist : float, 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. + + expandFactor : float, 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. + Raises ------ OSError @@ -160,6 +171,11 @@ def build_mapping_file(self, method='bilinear', # noqa: C901 raise ValueError(f'method {method} not supported for destination ' f'grid of type PointCollectionDescriptor.') + if (expandFactor is not None or expandDist is not None) and \ + method != 'conserve': + raise ValueError('expandFactor and expandDist only have an impact ' + 'with method=conserve') + if self.mappingFileName is None or \ os.path.exists(self.mappingFileName): # a valid weight file already exists, so nothing to do @@ -191,7 +207,9 @@ def build_mapping_file(self, method='bilinear', # noqa: C901 self.sourceDescriptor.to_scrip(sourceFileName) - self.destinationDescriptor.to_scrip(destinationFileName) + self.destinationDescriptor.to_scrip(destinationFileName, + expandDist=expandDist, + expandFactor=expandFactor) args = [rwgPath, '--source', sourceFileName,