From 36f3550e2f8ece589c6f5501202379e7e1597a66 Mon Sep 17 00:00:00 2001 From: Bin Xiao Date: Mon, 13 Jan 2020 22:07:00 +0800 Subject: [PATCH] Add initzonal_linear option in the initialization part which is alternative to the original initzonal option. The initzonal_linear option allow one to initialize the missing data with linearly-interpolated value along zonal direction, by which the converge iteration will be significantly reduced and the filled results will be more smooth in some circumstances. In addition, an initial_value can by used for the case in which most of the array are missing, and neither initzonal_linear nor initzonal option still work. --- gridfill/_gridfill.pyx | 106 +++++++++++++++++++++++++++++++++++------ gridfill/gridfill.py | 13 +++-- 2 files changed, 102 insertions(+), 17 deletions(-) diff --git a/gridfill/_gridfill.pyx b/gridfill/_gridfill.pyx index 415d273..48c9441 100644 --- a/gridfill/_gridfill.pyx +++ b/gridfill/_gridfill.pyx @@ -116,7 +116,10 @@ cdef void initialize_missing(double[:, :] grid, int[:, :] mask, unsigned int nlat, unsigned int nlon, - int initialize_zonal) nogil: + int initialize_zonal, + int initialize_zonal_linear, + int cyclic, + double initial_value) : """ Initialize the missing values in a grid in-place. @@ -143,11 +146,23 @@ cdef void initialize_missing(double[:, :] grid, values. If zero, use the value `0` as the initial guess for all missing values. + * initialize_zonal_linear [int] + If non-zero, take the zonal mean as the initial guess for missing + values. If zero, use the value `0` as the initial guess for all + missing values. """ - cdef unsigned int i, j, n - cdef double zonal_mean, initial_value - for i in range(nlat): - if initialize_zonal: + cdef unsigned int i, j, n, jm1, jp1, jj, njj, n_segment_s, n_segment_e + cdef double zonal_mean + cdef np.ndarray[UINT32_t, ndim=1] segment_s = np.zeros([nlon], + dtype=np.uint32)# x index of mask segment starters + cdef np.ndarray[FLOAT64_t, ndim=1] segment_s_val = np.zeros([nlon], + dtype=np.float64)# values of mask segment starters + cdef np.ndarray[UINT32_t, ndim=1] segment_e = np.zeros([nlon], + dtype=np.uint32)# x index of mask segment endings + cdef np.ndarray[FLOAT64_t, ndim=1] segment_e_val = np.zeros([nlon], + dtype=np.float64)# values of mask segment endings + if initialize_zonal: + for i in range(nlat): n = 0 zonal_mean = 0 for j in range(nlon): @@ -157,11 +172,58 @@ cdef void initialize_missing(double[:, :] grid, if n > 0: zonal_mean /= n initial_value = zonal_mean - else: - initial_value = 0 - for j in range(nlon): - if mask[i, j]: - grid[i, j] = initial_value + for j in range(nlon): + if mask[i, j]: + grid[i, j] = initial_value + elif initialize_zonal_linear: + if not cyclic : + raise ValueError('data must be cyclic when initialize_zonal_linear') + for i in range(nlat): + n_segment_s=0 + n_segment_e=0 + segment_s_val[:]=0. + segment_e_val[:]=0. + for j in range(nlon): + longitude_indices(j, nlon, cyclic, &jm1, &jp1) + if not mask[i,jm1] and mask[i,j]: #segment starters + segment_s[n_segment_s]=j + segment_s_val[n_segment_s] = grid[i,jm1] + n_segment_s += 1 + elif mask[i,jm1] and not mask[i,j]: #segment endings + segment_e[n_segment_e]=j + segment_e_val[n_segment_e] = grid[i,j] + n_segment_e += 1 + if n_segment_s != n_segment_e: + raise ValueError('n_segment_s not eq n_segment_e') + + if n_segment_s == 0: + for j in range(nlon): + if mask[i, j]: + grid[i, j] = initial_value + elif segment_s[0] < segment_e[0]: + for j in range(n_segment_s): + njj=0 + for jj in range(segment_s[j],segment_e[j]): + grid[i, jj] = njj * (segment_e_val[j] - segment_s_val[j])/(segment_e[j] - segment_s[j]) + segment_s_val[j] + njj += 1 + else: + for j in range(n_segment_s-1): + njj=0 + for jj in range(segment_s[j],segment_e[j+1]): + grid[i, jj] = njj * (segment_e_val[j+1] - segment_s_val[j])/(segment_e[j+1] - segment_s[j]) + segment_s_val[j] + njj += 1 + njj=0 + for jj in range(segment_s[n_segment_s-1],nlon): + grid[i, jj] = njj * (segment_e_val[0] - segment_s_val[n_segment_s-1])/(nlon-segment_s[n_segment_s-1] + segment_e[0]) + segment_s_val[n_segment_s-1] + njj += 1 + for jj in range(0,segment_e[0]): + grid[i, jj] = njj * (segment_e_val[0] - segment_s_val[n_segment_s-1])/(nlon-segment_s[n_segment_s-1] + segment_e[0]) + segment_s_val[n_segment_s-1] + njj += 1 + else: + for i in range(nlat): + for j in range(nlon): + if mask[i, j]: + grid[i, j] = initial_value cdef void poisson_fill(double[:, :] grid, @@ -173,8 +235,10 @@ cdef void poisson_fill(double[:, :] grid, unsigned int itermax, int cyclic, int initialize_zonal, + int initialize_zonal_linear, unsigned int *numiter, - double *resmax) nogil: + double *resmax, + double initial_value) : """ Fill missing values in a grid by iterative relaxation. @@ -215,6 +279,11 @@ cdef void poisson_fill(double[:, :] grid, values. If zero, use the value `0` as the initial guess for all missing values. + * initialize_zonal_linear [int] + If non-zero, take the zonal linear interpolation as the initial guess for missing + values. If zero, use the value `0` as the initial guess for all + missing values. + * numiter [unsigned int *] The number of iterations used to fill the grid. @@ -233,7 +302,7 @@ cdef void poisson_fill(double[:, :] grid, resmax[0] = 0 return # Set initial values for all missing values. - initialize_missing(grid, mask, nlat, nlon, initialize_zonal) + initialize_missing(grid, mask, nlat, nlon, initialize_zonal, initialize_zonal_linear, cyclic, initial_value) dp25 = 0.25 _numiter = 0 _resmax = 0 @@ -262,7 +331,9 @@ def poisson_fill_grids(double[:, :, :] grids, double tolerance, unsigned int itermax, int cyclic, - int initialize_zonal): + int initialize_zonal, + int initialize_zonal_linear, + double initial_value): """ Fill missing values in grids by iterative relaxation. @@ -297,6 +368,11 @@ def poisson_fill_grids(double[:, :, :] grids, values. If zero, use the value `0` as the initial guess for all missing values. + * initialize_zonal_linear [int] + If non-zero, take the zonal linear interpolation as the initial guess for missing + values. If zero, use the value `0` as the initial guess for all + missing values. + **Returns:** * numiter [numpy.ndarray[unsigned int, ndim=1] @@ -338,6 +414,8 @@ def poisson_fill_grids(double[:, :, :] grids, itermax, cyclic, initialize_zonal, + initialize_zonal_linear, &numiter[grid_num], - &resmax[grid_num]) + &resmax[grid_num], + initial_value) return (numiter, resmax) diff --git a/gridfill/gridfill.py b/gridfill/gridfill.py index ff2cfe3..a375b0c 100644 --- a/gridfill/gridfill.py +++ b/gridfill/gridfill.py @@ -64,8 +64,8 @@ def _recover_data(grid, info): return grid -def fill(grids, xdim, ydim, eps, relax=.6, itermax=100, initzonal=False, - cyclic=False, verbose=False): +def fill(grids, xdim, ydim, eps, relax=.6, itermax=100, initzonal=False, initzonal_linear=False, + cyclic=False, initial_value=0.0, verbose=False): """ Fill missing values in grids with values derived by solving Poisson's equation using a relaxation scheme. @@ -97,6 +97,11 @@ def fill(grids, xdim, ydim, eps, relax=.6, itermax=100, initzonal=False, missing values will be initialized to the zonal mean. Defaults to *False*. + *initzonal_linear* + If *False* missing values will be initialized to zero, if *True* + missing values will be initialized to the zonal mean. Defaults + to *False*. + *cyclic* Set to *False* if the x-coordinate of the grid is not cyclic, set to *True* if it is cyclic. Defaults to *False*. @@ -119,7 +124,9 @@ def fill(grids, xdim, ydim, eps, relax=.6, itermax=100, initzonal=False, # Call the computation subroutine: niter, resmax = _poisson_fill_grids(grids, masks, relax, eps, itermax, 1 if cyclic else 0, - 1 if initzonal else 0) + 1 if initzonal else 0, + 1 if initzonal_linear else 0, + initial_value) grids = _recover_data(grids, info) converged = np.logical_not(resmax > eps) # optional performance information: