Интерполяция трехмерных данных (неправильная вертикальная сетка) в регулярную вертикальную сетку. Циклы оптимизации - PullRequest
2 голосов
/ 08 марта 2012

Я работаю с 3D-полем.Горизонтальная сетка нерегулярна, а также вертикальная сетка нерегулярна и отличается от одного узла к другому в горизонтальной сетке (из-за сигма-вертикальных координат).

Я хочу интерполировать эти данные по регулярным вертикальным координатам, и сейчас единственный способ, который я нашел, - это перебрать каждую горизонтальную точку, чтобы затем выполнить простую одномерную вертикальную интерполяцию.Это легко реализовать (см. Ниже простой пример), но совершенно неэффективно (10 минут для интерполяции одного временного шага в моем исходном случае, и у меня более 10 000 временных шагов для последующей обработки) #!/ usr / bin / python импортировать numpy как np из scipy import интерполировать

## Horizontal mesh
x = range(0,300); y = range(0,200)
[X,Y] = np.meshgrid(x,y)

## Input 
z1   = np.array([0.,5.,10.,20.,50.,100.,200.,500.,1000.])
Z_in = np.zeros((z1.shape[0],X.shape[0],X.shape[1]))
for i in range(0,X.shape[0]):
  for j in range(0,X.shape[1]):
     Z_in[:,i,j]  = z1[:]
Z_in = Z_in+np.random.randint(-2,2,size=(z1.shape[0],X.shape[0],X.shape[1]))
var_in = np.random.randint(-10,10,size=(z1.shape[0],X.shape[0],X.shape[1]))

## Output
zout = np.array([0.,2.5,5.,7.5,10.,15.,20.,35.,50.,75.,100.,150.,200.,500.,1000.])
var_out = np.zeros((zout.shape[0],X.shape[0],X.shape[1]))

##Interpolation
for i in range(0,X.shape[0]):
   for j in range(0,X.shape[1]):
      var_out[:,i,j] = interpolate.interp1d(Z_in[:,i,j],var_in[:,i,j],bounds_error=False,fill_value=np.nan)(zout)

Если бы моя первоначальная вертикальная сетка была регулярной в каждой точке, это было бы легко в соответствии с предыдущим постом (http://stackoverflow.com/questions/7755871/interpolating-a-3d-array-in-python-how-to-avoid-for-loops), но в моемВ случае, если я не могу найти простое решение. Я попытался использовать сопоставление с функцией:

def interp_perso(Z_in,var_in,zout): return interpolate.interp1d(Z_in[:],var_in[:],bounds_error=False,fill_value=np.nan)(zout)

Zout =  np.zeros((X.shape[0],X.shape[1]))
for l in range(0,zout.shape[0]):
   Zout[:,:] = zout[l]
   var_out[l] = map(interp_perso,Z_in,var_in,Zout)

, но я все еще пришел к ошибке interp1d: "ValueError: массив x должен иметь ровно одно измерение."Вы знаете какой-либо другой способ реализовать это? Я что-то реализовал неправильно? Я был бы очень рад получить другие предложения.

Заранее спасибо.

...