Regrid 3DデカルトからポーラーへのNetCDFデータ

Question

私は時間、x、yの次元でNetCDF変数を持っています。現在はデカルト座標系ですが、極座標でデータが必要です。私はこれを行う関数を作成しようとしましたが、私はそれを正しくするように見えることはできません。これを行う簡単な方法があれば、誰にも分かりますか？現時点で

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は、サイズ176と軸1の境界の外にある」というエラーが表示されます。

何か助けていただければ幸いです。

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()