使用 Python Regrid 3D 笛卡尔到 Polar NetCDF 数据

Question

我有一个 NetCDF 变量，其时间维度为 x、y。它目前在笛卡尔坐标中，但我需要极坐标中的数据。我曾尝试创建一个函数来执行此操作，但我似乎无法做到这一点。有谁知道是否有更简单的方法可以做到这一点？

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,:,:])

目前，我收到一条错误消息，提示“索引 176 超出轴 1 的范围，大小为 176”，位于 * 行。

感谢任何帮助。

Answer 1

这是一种（未优化的）最近邻方法，用于将数据从常规笛卡尔网格重新映射到常规极坐标网格。坐标间距dr 和dphi 应根据您的需要进行调整。您还可以搜索 N 个最近的邻居并应用一些统计度量，例如距离加权平均值左右。对于大型数据集，我建议使用pykdtree 来加快最近邻搜索。对于重新映射地理空间数据，有一个名为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()