Можно ли векторизовать интерполяцию неравномерно разнесенного многомерного массива? - PullRequest
0 голосов
/ 23 марта 2019

В настоящее время у меня есть подпрограмма, которая берет 5D матрицу (x, y, z, t, p) и интерполирует данные в серию окружностей в пространстве x-y с разными радиусами. Для каждого p - p представляет модельное решение - полная интерполяция одного p (интерполяция x, y, z и t) занимает около 17 минут. К тому времени, когда я достигну p = 5, время ЦП будет превышено, а ядро ​​сброшено. Размер 5D матрицы составляет 1701,2201,40,48,6. Я не могу просто импортировать эти данные, поэтому у меня есть функция, которая берет первые 4 измерения и применяет операцию интерполяции. Я не могу использовать scipy.griddata, потому что интервал x-y (lat-lon) является неоднородным и изменяется в пространстве. Круги, на которые я пытаюсь интерполировать, имеют следующие радиусы (20, 30, 40, 50, 75, 100,125,150,175,200,250,300, 400,500,600,700) в километрах. Длина и ширина, простирающиеся от центра к краям окружностей, уже определены с помощью отдельного расчета. Каждый круг содержит около 72 точек с интервалом в 5 градусов от 0 до 360. Окончательный размер матрицы, если вычисление закончилось успешно, составит 6,48,17,72. Первые два индекса - это t и p, или время и результат модели, а 17 и 72 представляют номер круга, обозначенный радиусом, и точки вокруг круга. Моя рутина работает, если я делаю одно модельное решение, и я подтвердил, что результат реалистичен. Мне нужно найти способ ускорить этот расчет, чтобы избежать превышения времени процессора и сделать эту программу пригодной для будущего использования. Я просмотрел несколько вопросов, связанных с интерполяцией, и решил, что не могу использовать scipy.griddata; но мало информации, совпадающей с моим типом проблемы. Прямо сейчас я перебираю np.interp1d. Должно быть лучшее решение, возможно, что-то, связанное с векторизацией

У меня есть функция, которая берет файл данных (называемый File), импортирует его и использует сокращенную область ввода, а именно подмножество матрицы (например, lon-> lon [minrow: maxrow, mincol: maxcol]) ), чтобы сократить время вычислений при циклическом прохождении процедуры интерполяции. A

После того, как данные импортированы, я интерполирую через серию циклов (сетка неоднородна, и единственное, что, кажется, работает, это интерполяция измерения по измерению). Как видно из кода ниже, я перебираю каждое мгновение по времени (N), окружности (M) и угловому положению (P). Внутри Geointerp интерполяция дополнительно разбивается на серию циклов, основанных на высоте, длине и ширине прямоугольника. Внутри Geointerp я использую подпрограмму, которая идентифицирует самое близкое местоположение к интересующей точке, и я создаю квадратное подмножество (по сути, 3X3) и использую его в циклическом выражении lat и lon, а не по всему пространству lat-lon. Коды вставлены ниже ...

  1. Основная функция: отображение только соответствующих частей кода ....
N,M,P=Lonrange.shape #Lonrange is the lons of the circle (N=time,M=num of circles, P=num of angles)
n,m,p,q=var.shape #m=time index,n=height index,p=lon index, q=lat index
varinterp=np.zeros((N,M,P,m))
for j in range(N):
     for k in range(M):
          for l in range(P):
                  varinterp[j][k][l][:]=Geointerp(Lonrange[j][k][l],Latrange[j][k][l],Lon,Lat,var[j,:,:,:])
  1. Внутри Геоинтерп Примечание: помните, что Lon, Lat, Lonrange, Latrange и var являются входными данными
m,n=Lon.shape
lon=np.reshape(Lon,(m*n,1)) #changing the matrix to an array (speeds up calculation
lat=np.reshape(Lat,(m*n),1))
row,col=Minimum_Difference_2(Lon,Lat,lon,lat,Lonrange,Latrange)  #Finds row and col that is closest to the Lonrange and Latrange point!
LATs=Lat[row-1:row+2,col-1:col+2]
LONs=Lon[row-1:row+2,col-1:col+2]  #essentially creates a 3X3 subset
rr,cc=LONs.shape
M=var.shape
VAR=var[:,row-1:row+2,col-1:col+2] #
#Going through first dimesion, x
varinterp1=np.zeros((M[0],rr)) #var is a 3D matrix within Geointerp (z,x,y).  M[0]=the z index
lattest=np.zeros(rr)
for i in range(M[0]):
    for j in range(rr):
         if row==0 or col=0:
              varinterp1[i][j]=np.float('nan')
              lattest[j]=np.float('nan')
         else:
              varinterp1[i][j]=np.interp(Lonrange,LONs,VAR[i,j,:])
              lattest[j]=np.interp(Lonrange,LONs[j,:],LATs[j,:])
##Now next dimension (y)       
varinterp2=np.zeros((M[0]))
for i in range(M[0]):
        if row==0 or col==0: #If the row and col are 0, then outside interp range
               varinterp2[i]=np.float('nan')
        else:
               varinterp2[i]=np.interp(Latrange,lattest,varinterp1[i,:])

return varinterp2

ЭТО ДВУХ КУЗОВ КОДА, СОДЕРЖАЩИХ СООТВЕТСТВУЮЩУЮ ИНФОРМАЦИЮ. VARINTERP2 ДОЛЖЕН БЫТЬ ЗАКЛЮЧИТЕЛЬНЫМ РЕЗУЛЬТАТОМ, И ПРИ ОБЪЕДИНЕНИИ ДВУХ ВРЕМЕННЫХ ФУНКЦИЙ ДОЛЖЕН БЫТЬ РЕЗУЛЬТАТОВ В ЖЕЛАЕМ МАТРИЦЕ ВЫШЕ 6,48,17,72. Можно ли векторизовать какую-либо из процедур интерполяции, чтобы избежать чрезмерного зацикливания?

Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.
...