From 183cd8ce413c0c0d3c97f31bd1f4bbd24c20d7c9 Mon Sep 17 00:00:00 2001 From: Joseph Hamman Date: Fri, 22 Jun 2018 11:36:58 -0700 Subject: [PATCH 1/6] add mask to ds_in/ds_out --- xesmf/backend.py | 18 +++++++++++++++++- xesmf/frontend.py | 11 ++++++++++- xesmf/tests/test_frontend.py | 17 +++++++++++++++++ 3 files changed, 44 insertions(+), 2 deletions(-) diff --git a/xesmf/backend.py b/xesmf/backend.py index 2fa8d96..281df0f 100644 --- a/xesmf/backend.py +++ b/xesmf/backend.py @@ -52,7 +52,7 @@ def warn_lat_range(lat): warnings.warn("Latitude is outside of [-90, 90]") -def esmf_grid(lon, lat, periodic=False): +def esmf_grid(lon, lat, periodic=False, mask=None): ''' Create an ESMF.Grid object, for contrusting ESMF.Field and ESMF.Regrid @@ -70,6 +70,10 @@ def esmf_grid(lon, lat, periodic=False): Periodic in longitude? Default to False. Only useful for source grid. + mask : 2D numpy array, optional + Grid mask. Follows SCRIP convention where 1 is unmasked and 0 is + masked. + Returns ------- grid : ESMF.Grid object @@ -111,6 +115,18 @@ def esmf_grid(lon, lat, periodic=False): lon_pointer[...] = lon lat_pointer[...] = lat + # Follows SCRIP convention where 1 is unmasked and 0 is masked. + # See https://github.com/NCPP/ocgis/blob/61d88c60e9070215f28c1317221c2e074f8fb145/src/ocgis/regrid/base.py#L391-L404 + if mask is not None: + grid_mask = np.swapaxes(mask.astype(np.int32), 0, 1) + if not (grid_mask.shape == lon.shape): + raise ValueError("mask must have the same shape as the " + "latitude/longitude coordinates, got:" + "mask.shape = %s, lon.shape = %s" % (mask.shape, lon.shape)) + grid.add_item(ESMF.GridItem.MASK, staggerloc=ESMF.StaggerLoc.CENTER, + from_file=False) + grid.mask[0][:] = grid_mask + return grid diff --git a/xesmf/frontend.py b/xesmf/frontend.py index e2ced7d..5598adf 100644 --- a/xesmf/frontend.py +++ b/xesmf/frontend.py @@ -54,8 +54,14 @@ def ds_to_ESMFgrid(ds, need_bounds=False, periodic=None, append=None): lat = np.asarray(ds['lat']) lon, lat = as_2d_mesh(lon, lat) + if 'mask' in ds: + mask = np.asarray(ds['mask']) + print(mask.shape) + else: + mask = None + # tranpose the arrays so they become Fortran-ordered - grid = esmf_grid(lon.T, lat.T, periodic=periodic) + grid = esmf_grid(lon.T, lat.T, periodic=periodic, mask=mask) if need_bounds: lon_b = np.asarray(ds['lon_b']) @@ -83,6 +89,9 @@ def __init__(self, ds_in, ds_out, method, periodic=False, or 2D (Ny, Nx) for general curvilinear grids. Shape of bounds should be (N+1,) or (Ny+1, Nx+1). + If eiether dataset includes a 2d mask variable, that will also be + used to inform the regridding. + method : str Regridding method. Options are diff --git a/xesmf/tests/test_frontend.py b/xesmf/tests/test_frontend.py index 5b9d151..1609ebc 100644 --- a/xesmf/tests/test_frontend.py +++ b/xesmf/tests/test_frontend.py @@ -170,3 +170,20 @@ def test_regrid_with_1d_grid(): # clean-up regridder.clean_weight_file() + + +def test_build_regridder_with_masks(): + ds_in['mask'] = xr.DataArray( + np.random.randint(2, size=ds_in['data'].shape), + dims=('y', 'x')) + print(ds_in) + # 'patch' is too slow to test + for method in ['bilinear', 'conservative', 'nearest_s2d', 'nearest_d2s']: + regridder = xe.Regridder(ds_in, ds_out, method) + + # check screen output + assert repr(regridder) == str(regridder) + assert 'xESMF Regridder' in str(regridder) + assert method in str(regridder) + + regridder.clean_weight_file() From 0ba10d490a3de65ca35cc974f75c9d36961db5c5 Mon Sep 17 00:00:00 2001 From: Joseph Hamman Date: Fri, 22 Jun 2018 13:18:52 -0600 Subject: [PATCH 2/6] add mask values to Regrid call --- xesmf/backend.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/xesmf/backend.py b/xesmf/backend.py index 281df0f..ed3317a 100644 --- a/xesmf/backend.py +++ b/xesmf/backend.py @@ -256,7 +256,8 @@ def esmf_regrid_build(sourcegrid, destgrid, method, # if the destination grid is larger than the source grid. regrid = ESMF.Regrid(sourcefield, destfield, filename=filename, regrid_method=esmf_regrid_method, - unmapped_action=ESMF.UnmappedAction.IGNORE) + unmapped_action=ESMF.UnmappedAction.IGNORE, + src_mask_values=[0], dst_mask_values=[0]) return regrid From c39efe98af933ec5a288375284c64ad618b7f5ac Mon Sep 17 00:00:00 2001 From: Joseph Hamman Date: Fri, 22 Jun 2018 12:19:22 -0700 Subject: [PATCH 3/6] fix typo in doc string --- xesmf/frontend.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/xesmf/frontend.py b/xesmf/frontend.py index 5598adf..2b4565e 100644 --- a/xesmf/frontend.py +++ b/xesmf/frontend.py @@ -89,7 +89,7 @@ def __init__(self, ds_in, ds_out, method, periodic=False, or 2D (Ny, Nx) for general curvilinear grids. Shape of bounds should be (N+1,) or (Ny+1, Nx+1). - If eiether dataset includes a 2d mask variable, that will also be + If either dataset includes a 2d mask variable, that will also be used to inform the regridding. method : str From bc0d5e851bc2f46af59e0917251e8c0720f26256 Mon Sep 17 00:00:00 2001 From: Joseph Hamman Date: Fri, 22 Jun 2018 13:40:27 -0600 Subject: [PATCH 4/6] mask values are only 0 or 1 --- xesmf/backend.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/xesmf/backend.py b/xesmf/backend.py index ed3317a..df152e8 100644 --- a/xesmf/backend.py +++ b/xesmf/backend.py @@ -119,10 +119,12 @@ def esmf_grid(lon, lat, periodic=False, mask=None): # See https://github.com/NCPP/ocgis/blob/61d88c60e9070215f28c1317221c2e074f8fb145/src/ocgis/regrid/base.py#L391-L404 if mask is not None: grid_mask = np.swapaxes(mask.astype(np.int32), 0, 1) + grid_mask = np.where(grid_mask == 0, 0, 1) if not (grid_mask.shape == lon.shape): - raise ValueError("mask must have the same shape as the " - "latitude/longitude coordinates, got:" - "mask.shape = %s, lon.shape = %s" % (mask.shape, lon.shape)) + raise ValueError( + "mask must have the same shape as the latitude/longitude" + "coordinates, got: mask.shape = %s, lon.shape = %s" % + (mask.shape, lon.shape)) grid.add_item(ESMF.GridItem.MASK, staggerloc=ESMF.StaggerLoc.CENTER, from_file=False) grid.mask[0][:] = grid_mask From e2501c15dcd62db7a6798d0e2a72f9881037bba6 Mon Sep 17 00:00:00 2001 From: JiaweiZhuang Date: Mon, 6 Aug 2018 13:13:16 -0400 Subject: [PATCH 5/6] Add normalization option for conservative regridding See https://github.com/JiaweiZhuang/xESMF/issues/17. This together with masking can address the normalization issue in https://github.com/JiaweiZhuang/xESMF/pull/23 --- xesmf/backend.py | 10 ++++++++++ xesmf/frontend.py | 1 + 2 files changed, 11 insertions(+) diff --git a/xesmf/backend.py b/xesmf/backend.py index df152e8..b0d0b4a 100644 --- a/xesmf/backend.py +++ b/xesmf/backend.py @@ -193,6 +193,7 @@ def esmf_regrid_build(sourcegrid, destgrid, method, - 'bilinear' - 'conservative', **need grid corner information** + - 'conservative_normed', **need grid corner information** - 'patch' - 'nearest_s2d' - 'nearest_d2s' @@ -222,6 +223,7 @@ def esmf_regrid_build(sourcegrid, destgrid, method, # use shorter, clearer names for options in ESMF.RegridMethod method_dict = {'bilinear': ESMF.RegridMethod.BILINEAR, 'conservative': ESMF.RegridMethod.CONSERVE, + 'conservative_normed': ESMF.RegridMethod.CONSERVE, 'patch': ESMF.RegridMethod.PATCH, 'nearest_s2d': ESMF.RegridMethod.NEAREST_STOD, 'nearest_d2s': ESMF.RegridMethod.NEAREST_DTOS @@ -253,12 +255,20 @@ def esmf_regrid_build(sourcegrid, destgrid, method, assert not os.path.exists(filename), ( 'Weight file already exists! Please remove it or use a new name.') + # re-normalize conservative regridding results + # https://github.com/JiaweiZhuang/xESMF/issues/17 + if method == 'conservative_normed': + norm_type = ESMF.NormType.FRACAREA + else: + norm_type = ESMF.NormType.DSTAREA + # Calculate regridding weights. # Must set unmapped_action to IGNORE, otherwise the function will fail, # if the destination grid is larger than the source grid. regrid = ESMF.Regrid(sourcefield, destfield, filename=filename, regrid_method=esmf_regrid_method, unmapped_action=ESMF.UnmappedAction.IGNORE, + norm_type=norm_type, src_mask_values=[0], dst_mask_values=[0]) return regrid diff --git a/xesmf/frontend.py b/xesmf/frontend.py index 2b4565e..f4f5255 100644 --- a/xesmf/frontend.py +++ b/xesmf/frontend.py @@ -97,6 +97,7 @@ def __init__(self, ds_in, ds_out, method, periodic=False, - 'bilinear' - 'conservative', **need grid corner information** + - 'conservative_normed', **need grid corner information** - 'patch' - 'nearest_s2d' - 'nearest_d2s' From 3588a78c00cc3f9f0a9d313509c49fd7e3510db1 Mon Sep 17 00:00:00 2001 From: JiaweiZhuang Date: Mon, 6 Aug 2018 14:00:06 -0400 Subject: [PATCH 6/6] Pass corner infomation for 'conservative_normed' method --- xesmf/backend.py | 2 +- xesmf/frontend.py | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/xesmf/backend.py b/xesmf/backend.py index b0d0b4a..7e257bc 100644 --- a/xesmf/backend.py +++ b/xesmf/backend.py @@ -235,7 +235,7 @@ def esmf_regrid_build(sourcegrid, destgrid, method, '{}'.format(list(method_dict.keys()))) # conservative regridding needs cell corner information - if method == 'conservative': + if method in ['conservative', 'conservative_normed']: if not sourcegrid.has_corners: raise ValueError('source grid has no corner information. ' 'cannot use conservative regridding.') diff --git a/xesmf/frontend.py b/xesmf/frontend.py index f4f5255..aee87d0 100644 --- a/xesmf/frontend.py +++ b/xesmf/frontend.py @@ -83,7 +83,7 @@ def __init__(self, ds_in, ds_out, method, periodic=False, ds_in, ds_out : xarray DataSet, or dictionary Contain input and output grid coordinates. Look for variables ``lon``, ``lat``, and optionally ``lon_b``, ``lat_b`` for - conservative method. + add method. Shape can be 1D (Nlon,) and (Nlat,) for rectilinear grids, or 2D (Ny, Nx) for general curvilinear grids. @@ -125,7 +125,7 @@ def __init__(self, ds_in, ds_out, method, periodic=False, """ # record basic switches - if method == 'conservative': + if method in ['conservative', 'conservative_normed']: self.need_bounds = True periodic = False # bound shape will not be N+1 for periodic grid else: