CODEWITHSUNDEEP
pythonnumpyscipypython-xarray
alextc
alextc

Reputation: 3515

Use of scipy.interpolate.griddata for interpolation of data of multiple dimensions from xarray

I have a rainfall data as xarray DataSet named ds of three dimensions longitude, latitude and time:

<xarray.Dataset>
Dimensions:    (latitude: 691, longitude: 886, time: 1)
Coordinates:
  * longitude  (longitude) float64 112.0 112.0 112.1 112.1 ... 156.2 156.2 156.3
  * latitude   (latitude) float64 -9.975 -10.03 -10.08 ... -44.42 -44.47 -44.52
  * time       (time) datetime64[ns] 1980-01-03
Data variables:
    RAIN       (time, latitude, longitude) float64 0.0 0.0 0.0 ... 0.0 0.0 0.0

I would like to interpolate the rainfall values against another sets of longitudes and latitudes: EXAMPLE_FFDI_LON_XR_DA and EXAMPLE_FFDI_LAT_XR_DA. They have completely different values than the longitudes and latitudes of ds.

EXAMPLE_FFDI_LON_XR_DA:

<xarray.DataArray 'longitude' (longitude: 193)>
array([140.8    , 140.84792, ... ...], dtype=float32)
Coordinates:
  * longitude  (longitude) float32 140.8 140.84792 140.89584 ... 149.95209 150.0
Attributes:
    latIntersect:         0.0
    lonCentre:            145.4
    units:                degrees_east
    projectionType:       MERCATOR
    _CoordinateAxisType:  Lon

EXAMPLE_FFDI_LAT_XR_DA:

<xarray.DataArray 'latitude' (latitude: 106)>
array([-39.2     , -39.149525, ... ...], dtype=float32)
Coordinates:
  * latitude  (latitude) float32 -39.2 -39.149525 -39.09905 ... -33.950478 -33.9
Attributes:
    latIntersect:         0.0
    lonCentre:            145.4
    units:                degrees_north
    projectionType:       MERCATOR
    _CoordinateAxisType:  Lat

I thought of using xarray xarray.DataArray.interp function but this only supports the nearest method. I am new to scipy but thought it would better suit my need to interpolate using the scipy library scipy.interpolate.griddata function. How can I go about using this function with my data? A working example would be helpful.

Upvotes: 3

Views: 3150

Answers (1)

Marc
Marc

Reputation: 1629

I'm not sure it's possible to have scipy interact with the DataArray objects as you have them now. What you can do is the following:

from scipy.interpolate import griddata
import xarray
import numpy as np
# construct a dummy df like you have
values = np.random.rand(10, 30, 30) * 100
df = xarray.DataArray(values, [('time', range(10)), ('latitude', range(30)), ('longitude', range(30))])
# extract the values from the array
latitude = df['latitude'].values
longitude = df['longitude'].values
length_values = len(latitude) * len(longitude)

# create grid
x_grid, y_grid = np.meshgrid(latitude, longitude)
points = np.empty((length_values, 2))
points[:, 0] = x_grid.flatten()
points[:, 1] = y_grid.flatten()

# select a single timestep in which to interpolate the data
rain_at_locations = df.where(df.time == 3, drop=True)
values = np.array(rain_at_locations.values.flatten())

#your test locations
EXAMPLE_FFDI_LAT_XR_DA = range(5, 15) # extract these from you arrays
EXAMPLE_FFDI_LON_XR_DA = range(20, 30)
grid_x, grid_y = np.meshgrid(EXAMPLE_FFDI_LAT_XR_DA, EXAMPLE_FFDI_LON_XR_DA)

interpolated_values = griddata(points, values, (grid_x, grid_y), method='cubic')

These values can then be resupplied to a data array.

Upvotes: 2

Related Questions