Regrid 3D Cartesian에서 Polar로 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()