Andrew Dawson created a neat python module called gridfill. gridfill wraps a FORTRAN90 program that fills missing values in a grid by solving Poisson's equation using an iterative relaxation scheme. It is very similar o NCL's poisson_grid_fill.

Here is an example on how to use module.

In [3]:

import numpy.ma as ma
from scipy.io import loadmat

data = loadmat('./data/mercator_temperature.mat', squeeze_me=True)

x, y = data['x'], data['y']
bbox = [x.min(), x.max(), y.min(), y.max()]

temp = ma.masked_invalid(data['temp'][0, ...])

The data I loaded above are from an online example about Optimum Interpolation.

The next cell contains a simple function to create a map for the figures.

In [4]:

import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import cartopy.feature as cfeature
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER

LAND = cfeature.NaturalEarthFeature('physical', 'land', '50m',
                                    edgecolor='face',
                                    facecolor=cfeature.COLORS['land'])

def make_map(bbox, projection=ccrs.PlateCarree()):
    fig, ax = plt.subplots(figsize=(8, 6),
                           subplot_kw=dict(projection=projection))
    ax.set_extent(bbox)
    ax.add_feature(LAND, facecolor='0.75')
    ax.coastlines(resolution='50m')
    gl = ax.gridlines(draw_labels=True)
    gl.xlabels_top = gl.ylabels_right = False
    gl.xformatter = LONGITUDE_FORMATTER
    gl.yformatter = LATITUDE_FORMATTER
    return fig, ax

Here is the original data,

In [5]:

fig, ax = make_map(bbox=bbox)
cs = ax.pcolormesh(x, y, temp)
cbar = fig.colorbar(cs, extend='both', shrink=0.85)

There is nothing wrong with the data, but let's pretend it has several gaps due to clouds or sensor errors. (The functions below will help to create the gaps in the data.)

In [6]:

import numpy as np


def random_mask(data, frac=0.05):
    N = data.size
    frac = np.int(frac * N)
    # FIXME: Must excluded the land mask, so frac has a real meaning!
    idx = np.random.random_integers(0, N-1, frac)
    return idx

def mask_data(data, idx):
    out = data.copy()
    i, j = np.unravel_index(idx, data.shape)
    out[i, j] = np.NaN
    return ma.masked_invalid(out)

This is how our "incomplete" data look like.

In [7]:

idx = random_mask(temp, frac=0.05)
t05 = mask_data(temp, idx)

fig, ax = make_map(bbox)
ax.set_title('Removed {:0.2f} % of data.'.format(100 * idx.size / temp.size))
cs = ax.pcolormesh(x, y, t05)
cbar = fig.colorbar(cs, extend='both', shrink=0.85)

Now we can use gridfill to recover something that, hopefully, will look like the original data.

The arguments for fill are the gridded data, and the xdim and ydim dimensions that correspond to the x-coordinate and y-coordinate. The keyword arguments are:

eps: tolerance for the relaxation
relax: relaxation constant
itermax: maximum number of iterations
cyclic: True if thed ata is cyclic

In [8]:

from gridfill import fill

kw = dict(eps=1e-4, relax=0.6, itermax=1e4, initzonal=False,
          cyclic=False, verbose=True)

filled, converged = fill(t05, 1, 0, **kw)

if converged:
    fig, ax = make_map(bbox)
    cs = ax.pcolormesh(x, y, filled)
    cbar = fig.colorbar(cs, extend='both', shrink=0.85)

[0] relaxation converged (3745 iterations with maximum residual 9.991e-05)

Easy as pie! I used gridfill in some real applications with good results, like filling gaps in HF radar currents and to creating gridded CTD data for scientific cruises.

I am trying to expand gridfill so it does not fill over a certain mask, like continents or a bathymetric line. As soon as I have something that works I will create another post.

I almost forgot! gridfill has a fill_cube alternative that operates on iris cubes.

In [9]:

HTML(html)

Out[9]:

This post was written as an IPython notebook. It is available for download or as a static html.

python4oceanographers by Filipe Fernandes is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License.
Based on a work at https://ocefpaf.github.io/.

python4oceanographers

Turning ripples into waves

Filling gaps in gridded data

Comments