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 initzonal_linear option in the initialization part #19

Open
wants to merge 1 commit into
base: main
Choose a base branch
from
Open

Add initzonal_linear option in the initialization part #19

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
106 changes: 92 additions & 14 deletions gridfill/_gridfill.pyx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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.

Expand All @@ -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):
Expand All @@ -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,
Expand All @@ -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.

Expand Down Expand Up @@ -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.

Expand All @@ -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
Expand Down Expand Up @@ -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.

Expand Down Expand Up @@ -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]
Expand Down Expand Up @@ -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)
13 changes: 10 additions & 3 deletions gridfill/gridfill.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down Expand Up @@ -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*.
Expand All @@ -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:
Expand Down