Interpolating data

xarray offers flexible interpolation routines, which have a similar interface to our indexing.

Note

interp requires scipy installed.

Scalar and 1-dimensional interpolation

Interpolating a DataArray works mostly like labeled indexing of a DataArray,

In [1]: da = xr.DataArray(np.sin(0.3 * np.arange(12).reshape(4, 3)),
   ...:                   [('time', np.arange(4)),
   ...:                    ('space', [0.1, 0.2, 0.3])])
   ...: 

# label lookup
In [2]: da.sel(time=3)
Out[2]: 
<xarray.DataArray (space: 3)>
array([ 0.42738 ,  0.14112 , -0.157746])
Coordinates:
    time     int64 3
  * space    (space) float64 0.1 0.2 0.3

# interpolation
In [3]: da.interp(time=2.5)
Out[3]: 
<xarray.DataArray (space: 3)>
array([0.700614, 0.502165, 0.258859])
Coordinates:
  * space    (space) float64 0.1 0.2 0.3
    time     float64 2.5

Similar to the indexing, interp() also accepts an array-like, which gives the interpolated result as an array.

# label lookup
In [4]: da.sel(time=[2, 3])
Out[4]: 
<xarray.DataArray (time: 2, space: 3)>
array([[ 0.973848,  0.863209,  0.675463],
       [ 0.42738 ,  0.14112 , -0.157746]])
Coordinates:
  * time     (time) int64 2 3
  * space    (space) float64 0.1 0.2 0.3

# interpolation
In [5]: da.interp(time=[2.5, 3.5])
Out[5]: 
<xarray.DataArray (time: 2, space: 3)>
array([[0.700614, 0.502165, 0.258859],
       [     nan,      nan,      nan]])
Coordinates:
  * space    (space) float64 0.1 0.2 0.3
  * time     (time) float64 2.5 3.5

To interpolate data with a numpy.datetime64() coordinate you can pass a string.

In [6]: da_dt64 = xr.DataArray([1, 3],
   ...:                        [('time', pd.date_range('1/1/2000', '1/3/2000', periods=2))])
   ...: 

In [7]: da_dt64.interp(time='2000-01-02')
Out[7]: 
<xarray.DataArray ()>
array(2.)
Coordinates:
    time     datetime64[ns] 2000-01-02

The interpolated data can be merged into the original DataArray by specifying the time periods required.

In [8]: da_dt64.interp(time=pd.date_range('1/1/2000', '1/3/2000', periods=3))
Out[8]: 
<xarray.DataArray (time: 3)>
array([1., 2., 3.])
Coordinates:
  * time     (time) datetime64[ns] 2000-01-01 2000-01-02 2000-01-03

Interpolation of data indexed by a CFTimeIndex is also allowed. See Non-standard calendars and dates outside the Timestamp-valid range for examples.

Note

Currently, our interpolation only works for regular grids. Therefore, similarly to sel(), only 1D coordinates along a dimension can be used as the original coordinate to be interpolated.

Multi-dimensional Interpolation

Like sel(), interp() accepts multiple coordinates. In this case, multidimensional interpolation is carried out.

# label lookup
In [9]: da.sel(time=2, space=0.1)
Out[9]: 
<xarray.DataArray ()>
array(0.973848)
Coordinates:
    time     int64 2
    space    float64 0.1

# interpolation
In [10]: da.interp(time=2.5, space=0.15)
Out[10]: 
<xarray.DataArray ()>
array(0.601389)
Coordinates:
    time     float64 2.5
    space    float64 0.15

Array-like coordinates are also accepted:

# label lookup
In [11]: da.sel(time=[2, 3], space=[0.1, 0.2])
Out[11]: 
<xarray.DataArray (time: 2, space: 2)>
array([[0.973848, 0.863209],
       [0.42738 , 0.14112 ]])
Coordinates:
  * time     (time) int64 2 3
  * space    (space) float64 0.1 0.2

# interpolation
In [12]: da.interp(time=[1.5, 2.5], space=[0.15, 0.25])
Out[12]: 
<xarray.DataArray (time: 2, space: 2)>
array([[0.888106, 0.867052],
       [0.601389, 0.380512]])
Coordinates:
  * time     (time) float64 1.5 2.5
  * space    (space) float64 0.15 0.25

interp_like() method is a useful shortcut. This method interpolates an xarray object onto the coordinates of another xarray object. For example, if we want to compute the difference between two DataArray s (da and other) staying on slightly different coordinates,

In [13]: other = xr.DataArray(np.sin(0.4 * np.arange(9).reshape(3, 3)),
   ....:                      [('time', [0.9, 1.9, 2.9]),
   ....:                      ('space', [0.15, 0.25, 0.35])])
   ....: 

it might be a good idea to first interpolate da so that it will stay on the same coordinates of other, and then subtract it. interp_like() can be used for such a case,

# interpolate da along other's coordinates
In [14]: interpolated = da.interp_like(other)

In [15]: interpolated
Out[15]: 
<xarray.DataArray (time: 3, space: 3)>
array([[0.786691, 0.911298,      nan],
       [0.912444, 0.788879,      nan],
       [0.347678, 0.069452,      nan]])
Coordinates:
  * time     (time) float64 0.9 1.9 2.9
  * space    (space) float64 0.15 0.25 0.35

It is now possible to safely compute the difference other - interpolated.

Interpolation methods

We use scipy.interpolate.interp1d() for 1-dimensional interpolation and scipy.interpolate.interpn() for multi-dimensional interpolation.

The interpolation method can be specified by the optional method argument.

In [16]: da = xr.DataArray(np.sin(np.linspace(0, 2 * np.pi, 10)), dims='x',
   ....:                   coords={'x': np.linspace(0, 1, 10)})
   ....: 

In [17]: da.plot.line('o', label='original')
Out[17]: [<matplotlib.lines.Line2D at 0x74fd5858bf98>]

In [18]: da.interp(x=np.linspace(0, 1, 100)).plot.line(label='linear (default)')
Out[18]: [<matplotlib.lines.Line2D at 0x74fd585c5d30>]

In [19]: da.interp(x=np.linspace(0, 1, 100), method='cubic').plot.line(label='cubic')
Out[19]: [<matplotlib.lines.Line2D at 0x74fd585c5eb8>]

In [20]: plt.legend()
Out[20]: <matplotlib.legend.Legend at 0x74fd5840d630>
_images/interpolation_sample1.png

Additional keyword arguments can be passed to scipy’s functions.

# fill 0 for the outside of the original coordinates.
In [21]: da.interp(x=np.linspace(-0.5, 1.5, 10), kwargs={'fill_value': 0.0})
Out[21]: 
<xarray.DataArray (x: 10)>
array([ 0.      ,  0.      ,  0.      ,  0.813798,  0.604023, -0.604023,
       -0.813798,  0.      ,  0.      ,  0.      ])
Coordinates:
  * x        (x) float64 -0.5 -0.2778 -0.05556 0.1667 ... 0.8333 1.056 1.278 1.5

# extrapolation
In [22]: da.interp(x=np.linspace(-0.5, 1.5, 10), kwargs={'fill_value': 'extrapolate'})
Out[22]: 
<xarray.DataArray (x: 10)>
array([-2.892544, -1.606969, -0.321394,  0.813798,  0.604023, -0.604023,
       -0.813798,  0.321394,  1.606969,  2.892544])
Coordinates:
  * x        (x) float64 -0.5 -0.2778 -0.05556 0.1667 ... 0.8333 1.056 1.278 1.5

Advanced Interpolation

interp() accepts DataArray as similar to sel(), which enables us more advanced interpolation. Based on the dimension of the new coordinate passed to interp(), the dimension of the result are determined.

For example, if you want to interpolate a two dimensional array along a particular dimension, as illustrated below, you can pass two 1-dimensional DataArray s with a common dimension as new coordinate.

advanced indexing and interpolation

For example:

In [23]: da = xr.DataArray(np.sin(0.3 * np.arange(20).reshape(5, 4)),
   ....:                   [('x', np.arange(5)),
   ....:                    ('y', [0.1, 0.2, 0.3, 0.4])])
   ....: 

# advanced indexing
In [24]: x = xr.DataArray([0, 2, 4], dims='z')

In [25]: y = xr.DataArray([0.1, 0.2, 0.3], dims='z')

In [26]: da.sel(x=x, y=y)
Out[26]: 
<xarray.DataArray (z: 3)>
array([ 0.      ,  0.42738 , -0.772764])
Coordinates:
    x        (z) int64 0 2 4
    y        (z) float64 0.1 0.2 0.3
Dimensions without coordinates: z

# advanced interpolation
In [27]: x = xr.DataArray([0.5, 1.5, 2.5], dims='z')

In [28]: y = xr.DataArray([0.15, 0.25, 0.35], dims='z')

In [29]: da.interp(x=x, y=y)
Out[29]: 
<xarray.DataArray (z: 3)>
array([ 0.556264,  0.634961, -0.466433])
Coordinates:
    x        (z) float64 0.5 1.5 2.5
    y        (z) float64 0.15 0.25 0.35
Dimensions without coordinates: z

where values on the original coordinates (x, y) = ((0.5, 0.15), (1.5, 0.25), (2.5, 0.35)) are obtained by the 2-dimensional interpolation and mapped along a new dimension z.

If you want to add a coordinate to the new dimension z, you can supply DataArray s with a coordinate,

In [30]: x = xr.DataArray([0.5, 1.5, 2.5], dims='z', coords={'z': ['a', 'b','c']})

In [31]: y = xr.DataArray([0.15, 0.25, 0.35], dims='z',
   ....:                  coords={'z': ['a', 'b','c']})
   ....: 

In [32]: da.interp(x=x, y=y)
Out[32]: 
<xarray.DataArray (z: 3)>
array([ 0.556264,  0.634961, -0.466433])
Coordinates:
    x        (z) float64 0.5 1.5 2.5
    y        (z) float64 0.15 0.25 0.35
  * z        (z) <U1 'a' 'b' 'c'

For the details of the advanced indexing, see more advanced indexing.

Interpolating arrays with NaN

Our interp() works with arrays with NaN the same way that scipy.interpolate.interp1d and scipy.interpolate.interpn do. linear and nearest methods return arrays including NaN, while other methods such as cubic or quadratic return all NaN arrays.

In [33]: da = xr.DataArray([0, 2, np.nan, 3, 3.25], dims='x',
   ....:                   coords={'x': range(5)})
   ....: 

In [34]: da.interp(x=[0.5, 1.5, 2.5])
Out[34]: 
<xarray.DataArray (x: 3)>
array([ 1., nan, nan])
Coordinates:
  * x        (x) float64 0.5 1.5 2.5

In [35]: da.interp(x=[0.5, 1.5, 2.5], method='cubic')
Out[35]: 
<xarray.DataArray (x: 3)>
array([nan, nan, nan])
Coordinates:
  * x        (x) float64 0.5 1.5 2.5

To avoid this, you can drop NaN by dropna(), and then make the interpolation

In [36]: dropped = da.dropna('x')

In [37]: dropped
Out[37]: 
<xarray.DataArray (x: 4)>
array([0.  , 2.  , 3.  , 3.25])
Coordinates:
  * x        (x) int64 0 1 3 4

In [38]: dropped.interp(x=[0.5, 1.5, 2.5], method='cubic')
Out[38]: 
<xarray.DataArray (x: 3)>
array([1.190104, 2.507812, 2.929688])
Coordinates:
  * x        (x) float64 0.5 1.5 2.5

If NaNs are distributed randomly in your multidimensional array, dropping all the columns containing more than one NaNs by dropna() may lose a significant amount of information. In such a case, you can fill NaN by interpolate_na(), which is similar to pandas.Series.interpolate().

In [39]: filled = da.interpolate_na(dim='x')

In [40]: filled
Out[40]: 
<xarray.DataArray (x: 5)>
array([0.  , 2.  , 2.5 , 3.  , 3.25])
Coordinates:
  * x        (x) int64 0 1 2 3 4

This fills NaN by interpolating along the specified dimension. After filling NaNs, you can interpolate:

In [41]: filled.interp(x=[0.5, 1.5, 2.5], method='cubic')
Out[41]: 
<xarray.DataArray (x: 3)>
array([1.308594, 2.316406, 2.738281])
Coordinates:
  * x        (x) float64 0.5 1.5 2.5

For the details of interpolate_na(), see Missing values.

Example

Let’s see how interp() works on real data.

# Raw data
In [42]: ds = xr.tutorial.open_dataset('air_temperature').isel(time=0)
---------------------------------------------------------------------------
gaierror                                  Traceback (most recent call last)
/usr/lib/python3.7/urllib/request.py in do_open(self, http_class, req, **http_conn_args)
   1323                 h.request(req.get_method(), req.selector, req.data, headers,
-> 1324                           encode_chunked=req.has_header('Transfer-encoding'))
   1325             except OSError as err: # timeout error

/usr/lib/python3.7/http/client.py in request(self, method, url, body, headers, encode_chunked)
   1243         """Send a complete request to the server."""
-> 1244         self._send_request(method, url, body, headers, encode_chunked)
   1245 

/usr/lib/python3.7/http/client.py in _send_request(self, method, url, body, headers, encode_chunked)
   1289             body = _encode(body, 'body')
-> 1290         self.endheaders(body, encode_chunked=encode_chunked)
   1291 

/usr/lib/python3.7/http/client.py in endheaders(self, message_body, encode_chunked)
   1238             raise CannotSendHeader()
-> 1239         self._send_output(message_body, encode_chunked=encode_chunked)
   1240 

/usr/lib/python3.7/http/client.py in _send_output(self, message_body, encode_chunked)
   1025         del self._buffer[:]
-> 1026         self.send(msg)
   1027 

/usr/lib/python3.7/http/client.py in send(self, data)
    965             if self.auto_open:
--> 966                 self.connect()
    967             else:

/usr/lib/python3.7/http/client.py in connect(self)
   1398 
-> 1399             super().connect()
   1400 

/usr/lib/python3.7/http/client.py in connect(self)
    937         self.sock = self._create_connection(
--> 938             (self.host,self.port), self.timeout, self.source_address)
    939         self.sock.setsockopt(socket.IPPROTO_TCP, socket.TCP_NODELAY, 1)

/usr/lib/python3.7/socket.py in create_connection(address, timeout, source_address)
    706     err = None
--> 707     for res in getaddrinfo(host, port, 0, SOCK_STREAM):
    708         af, socktype, proto, canonname, sa = res

/usr/lib/python3.7/socket.py in getaddrinfo(host, port, family, type, proto, flags)
    747     addrlist = []
--> 748     for res in _socket.getaddrinfo(host, port, family, type, proto, flags):
    749         af, socktype, proto, canonname, sa = res

gaierror: [Errno -3] Temporary failure in name resolution

During handling of the above exception, another exception occurred:

URLError                                  Traceback (most recent call last)
<ipython-input-42-13a3a9c99ae5> in <module>()
----> 1 ds = xr.tutorial.open_dataset('air_temperature').isel(time=0)

/opt/build/python-xarray-0.11.3/xarray/tutorial.py in open_dataset(name, cache, cache_dir, github_url, branch, **kws)
     69 
     70         url = '/'.join((github_url, 'raw', branch, fullname))
---> 71         _urlretrieve(url, localfile)
     72         url = '/'.join((github_url, 'raw', branch, md5name))
     73         _urlretrieve(url, md5file)

/usr/lib/python3.7/urllib/request.py in urlretrieve(url, filename, reporthook, data)
    245     url_type, path = splittype(url)
    246 
--> 247     with contextlib.closing(urlopen(url, data)) as fp:
    248         headers = fp.info()
    249 

/usr/lib/python3.7/urllib/request.py in urlopen(url, data, timeout, cafile, capath, cadefault, context)
    220     else:
    221         opener = _opener
--> 222     return opener.open(url, data, timeout)
    223 
    224 def install_opener(opener):

/usr/lib/python3.7/urllib/request.py in open(self, fullurl, data, timeout)
    523             req = meth(req)
    524 
--> 525         response = self._open(req, data)
    526 
    527         # post-process response

/usr/lib/python3.7/urllib/request.py in _open(self, req, data)
    541         protocol = req.type
    542         result = self._call_chain(self.handle_open, protocol, protocol +
--> 543                                   '_open', req)
    544         if result:
    545             return result

/usr/lib/python3.7/urllib/request.py in _call_chain(self, chain, kind, meth_name, *args)
    501         for handler in handlers:
    502             func = getattr(handler, meth_name)
--> 503             result = func(*args)
    504             if result is not None:
    505                 return result

/usr/lib/python3.7/urllib/request.py in https_open(self, req)
   1365         def https_open(self, req):
   1366             return self.do_open(http.client.HTTPSConnection, req,
-> 1367                 context=self._context, check_hostname=self._check_hostname)
   1368 
   1369         https_request = AbstractHTTPHandler.do_request_

/usr/lib/python3.7/urllib/request.py in do_open(self, http_class, req, **http_conn_args)
   1324                           encode_chunked=req.has_header('Transfer-encoding'))
   1325             except OSError as err: # timeout error
-> 1326                 raise URLError(err)
   1327             r = h.getresponse()
   1328         except:

URLError: <urlopen error [Errno -3] Temporary failure in name resolution>

In [43]: fig, axes = plt.subplots(ncols=2, figsize=(10, 4))

In [44]: ds.air.plot(ax=axes[0])
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-44-cb8f083667be> in <module>()
----> 1 ds.air.plot(ax=axes[0])

/opt/build/python-xarray-0.11.3/xarray/core/common.py in __getattr__(self, name)
    177                     return source[name]
    178         raise AttributeError("%r object has no attribute %r" %
--> 179                              (type(self).__name__, name))
    180 
    181     def __setattr__(self, name, value):

AttributeError: 'Dataset' object has no attribute 'air'

In [45]: axes[0].set_title('Raw data')
Out[45]: Text(0.5, 1.0, 'Raw data')

# Interpolated data
In [46]: new_lon = np.linspace(ds.lon[0], ds.lon[-1], ds.dims['lon'] * 4)
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-46-8ba5d53673bb> in <module>()
----> 1 new_lon = np.linspace(ds.lon[0], ds.lon[-1], ds.dims['lon'] * 4)

/opt/build/python-xarray-0.11.3/xarray/core/common.py in __getattr__(self, name)
    177                     return source[name]
    178         raise AttributeError("%r object has no attribute %r" %
--> 179                              (type(self).__name__, name))
    180 
    181     def __setattr__(self, name, value):

AttributeError: 'Dataset' object has no attribute 'lon'

In [47]: new_lat = np.linspace(ds.lat[0], ds.lat[-1], ds.dims['lat'] * 4)
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-47-5da4462467cd> in <module>()
----> 1 new_lat = np.linspace(ds.lat[0], ds.lat[-1], ds.dims['lat'] * 4)

/opt/build/python-xarray-0.11.3/xarray/core/common.py in __getattr__(self, name)
    177                     return source[name]
    178         raise AttributeError("%r object has no attribute %r" %
--> 179                              (type(self).__name__, name))
    180 
    181     def __setattr__(self, name, value):

AttributeError: 'Dataset' object has no attribute 'lat'

In [48]: dsi = ds.interp(lat=new_lat, lon=new_lon)
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-48-589e3eedfe36> in <module>()
----> 1 dsi = ds.interp(lat=new_lat, lon=new_lon)

NameError: name 'new_lat' is not defined

In [49]: dsi.air.plot(ax=axes[1])
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-49-433e5c8b0d13> in <module>()
----> 1 dsi.air.plot(ax=axes[1])

NameError: name 'dsi' is not defined

In [50]: axes[1].set_title('Interpolated data')
Out[50]: Text(0.5, 1.0, 'Interpolated data')
_images/interpolation_sample3.png

Our advanced interpolation can be used to remap the data to the new coordinate. Consider the new coordinates x and z on the two dimensional plane. The remapping can be done as follows

# new coordinate
In [51]: x = np.linspace(240, 300, 100)

In [52]: z = np.linspace(20, 70, 100)

# relation between new and original coordinates
In [53]: lat = xr.DataArray(z, dims=['z'], coords={'z': z})

In [54]: lon = xr.DataArray((x[:, np.newaxis]-270)/np.cos(z*np.pi/180)+270,
   ....:                    dims=['x', 'z'], coords={'x': x, 'z': z})
   ....: 

In [55]: fig, axes = plt.subplots(ncols=2, figsize=(10, 4))

In [56]: ds.air.plot(ax=axes[0])
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-56-cb8f083667be> in <module>()
----> 1 ds.air.plot(ax=axes[0])

/opt/build/python-xarray-0.11.3/xarray/core/common.py in __getattr__(self, name)
    177                     return source[name]
    178         raise AttributeError("%r object has no attribute %r" %
--> 179                              (type(self).__name__, name))
    180 
    181     def __setattr__(self, name, value):

AttributeError: 'Dataset' object has no attribute 'air'

# draw the new coordinate on the original coordinates.
In [57]: for idx in [0, 33, 66, 99]:
   ....:     axes[0].plot(lon.isel(x=idx), lat, '--k')
   ....: 

In [58]: for idx in [0, 33, 66, 99]:
   ....:     axes[0].plot(*xr.broadcast(lon.isel(z=idx), lat.isel(z=idx)), '--k')
   ....: 

In [59]: axes[0].set_title('Raw data')
Out[59]: Text(0.5, 1.0, 'Raw data')

In [60]: dsi = ds.interp(lon=lon, lat=lat)
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-60-eb48b9af7203> in <module>()
----> 1 dsi = ds.interp(lon=lon, lat=lat)

/opt/build/python-xarray-0.11.3/xarray/core/dataset.py in interp(self, coords, method, assume_sorted, kwargs, **coords_kwargs)
   1959 
   1960         coords = either_dict_or_kwargs(coords, coords_kwargs, 'interp')
-> 1961         indexers = OrderedDict(self._validate_indexers(coords))
   1962 
   1963         obj = self if assume_sorted else self.sortby([k for k in coords])

/opt/build/python-xarray-0.11.3/xarray/core/dataset.py in _validate_indexers(self, indexers)
   1407         invalid = [k for k in indexers if k not in self.dims]
   1408         if invalid:
-> 1409             raise ValueError("dimensions %r do not exist" % invalid)
   1410 
   1411         # all indexers should be int, slice, np.ndarrays, or Variable

ValueError: dimensions ['lon', 'lat'] do not exist

In [61]: dsi.air.plot(ax=axes[1])
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-61-433e5c8b0d13> in <module>()
----> 1 dsi.air.plot(ax=axes[1])

NameError: name 'dsi' is not defined

In [62]: axes[1].set_title('Remapped data')
Out[62]: Text(0.5, 1.0, 'Remapped data')
_images/interpolation_sample4.png