Regrid 3D Cartesian to Polar NetCDF Data with Python

Question

I have a NetCDF variable with the dimensions of time, x, y. It is currently in cartesian coordinates but I need the data in polar coordinates. I have tried creating a function to do this, but I cannot seem to get it right. Does anyone know if there is a simpler way to do this?

def regrid(x,y,xcent,ycent,vardat):
    x=np.subtract(x,xcent)
    y=np.subtract(y,ycent)
    threshmin = np.min(vardat)
    threshmax = np.max(vardat)
    rmax = np.ceil(np.sqrt(((x[-1]-x[0])/2.)**2 + ((y[-1]-y[0])/2.)**2))
    r = np.arange(0,rmax,(x[1]-x[0]))
    theta_inc = np.floor(np.arctan2(y[1]-y[0],(x[-1]-x[0])/2.)/np.pi*180.)
    if theta_inc <1.0:
        theta_inc = 1
    theta = np.arange(0,(360-theta_inc),theta_inc)
    r2d, theta2d = np.meshgrid(r,theta)
    x_polar = r2d*np.cos(np.pi/180.*theta2d)
    y_polar = r2d*np.sin(np.pi/180.*theta2d)
    x_range = np.arange(x[0],x[-1]+1,(x[1]-x[0]))
    y_range = np.arange(y[0],y[-1]+1,(y[1]-y[0]))
    field_rt = np.zeros((len(r),len(theta)))


    field_interp = interp2d(x_range,y_range,vardat,kind='linear')
    for i in np.arange(0,len(r)):
        for j in np.arange(0,len(theta)):

    *       field_rt[i,j] = field_interp(x_polar[i,j],y_polar[i,j])

    return r, theta, field_rt

r1,theta1, field = regrid(we_ea,no_so,124,124,olr[0,:,:])

At the moment, I get an error saying "index 176 is out of bounds for axis 1 with size 176" at the line with the *.

Any help is appreciated.

Answer 1

Here is a (not optimized) nearest neighbour approach to remap data from a regular cartesian grid to a regular polar grid. The coordinate spacings dr and dphi should be adjusted to your needs. You might also search for the N nearest neighbours and apply some statistical measure, e.g. a distance weighted average or so. For large datasets I would recommend using pykdtree to speed up the nearest neighbour search. For remapping geospatial data there is a nice and fast package called pyresample.

import numpy as np
import matplotlib.pyplot as plt

deg2rad = np.pi/180.0

# Define properties of cartesian grid
x_vals = np.arange(-1, 1, 0.01)
y_vals = np.arange(-1, 1, 0.01)
mx, my = np.meshgrid(y_vals, x_vals)

# Define data on cartesian grid
data_cart = np.sin(mx) + np.cos(my)

# Define properties of polar grid
dr = 0.1
dphi = 1*deg2rad
rmax = np.sqrt(x_vals.max()**2 + y_vals.max()**2)
r_vals = np.arange(0, rmax, dr)
phi_vals = np.arange(0, 2*np.pi, dphi)
if len(r_vals)*len(phi_vals) > len(x_vals)*len(y_vals):
    print "Warning: Oversampling"
mr, mphi = np.meshgrid(phi_vals, r_vals)

# Initialize data on polar grid with fill values
fill_value = -9999.0
data_polar = fill_value*np.ones((len(r_vals), len(phi_vals)))

# Define radius of influence. A nearest neighbour outside this radius will not
# be taken into account.
radius_of_influence = np.sqrt(0.1**2 + 0.1**2)

# For each cell in the polar grid, find the nearest neighbour in the cartesian
# grid. If it lies within the radius of influence, transfer the corresponding
# data.
for r, row_polar in zip(r_vals, range(len(r_vals))):
    for phi, col_polar in zip(phi_vals, range(len(phi_vals))):
        # Transform polar to cartesian
        x = r*np.cos(phi)
        y = r*np.sin(phi)

        # Find nearest neighbour in cartesian grid
        d = np.sqrt((x-mx)**2 + (y-my)**2)
        nn_row_cart, nn_col_cart = np.unravel_index(np.argmin(d), d.shape)
        dmin = d[nn_row_cart, nn_col_cart]

        # Transfer data
        if dmin <= radius_of_influence:
            data_polar[row_polar, col_polar] = data_cart[nn_row_cart, nn_col_cart]

# Mask remaining fill values
data_polar = np.ma.masked_equal(data_polar, fill_value)

# Plot results
plt.figure()
im = plt.pcolormesh(mx, my, data_cart)
plt.xlabel('x')
plt.ylabel('y')
plt.title('Cartesian')
plt.colorbar(im)

plt.figure()
ax = plt.subplot(111, projection='polar')
im = ax.pcolormesh(mr, mphi, data_polar)
ax.set_title('Polar')
plt.colorbar(im)

plt.show()