Существуют ли методы, которые выполняют так же, как использование "столбца стека" в Numpy PYTHON

Question

Я использую Python 2.7. Моя система работает под управлением Windows Vista, 32 бита.

У меня есть фрагмент кода, который читает яркость, широту и долготу и файл изображения (в расширении hdf). А затем попробуйте выполнить примерный ближайший сосед и отобразить его. Но когда он пытается приблизиться к ближайшему соседу, он выдает ошибку памяти.

Размер одного hdf-файла составляет 4,70 МБ, что, похоже, не слишком велико.

Вот мой код:

if __name__=="__main__":

    filename = ... ( the hdf file I have)      
    cumData, z = readAIRS_L1_VIS(filename)

    x, y = get_lat_lon(filename)  

    x0, xn = int(x.min()+1), int(x.max())
    y0, yn = int(y.min()+1), int(y.max())

    ncol = xn - x0 + 1
    nrow = yn - y0 + 1

    X, Y = np.meshgrid(np.arange(x0, xn+1), np.arange(y0, yn+1))
    img = interp_knn(np.column_stack((x.ravel(), y.ravel())),
            z.ravel(), np.column_stack((X.ravel(), Y.ravel())))
    img.shape = (nrow, ncol)

Тогда мои функции и импорт:

from pyhdf.SD import SD
import scipy as sc
import numpy as np
import pylab, os
import pyproj as proj
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
import scikits.ann as ann

def readAIRS_L1_VIS(filename,variable=None):
    allz=[]
    """
    function
        read hdf file for AIR Level 1B VIS
    input : AIRS HDF file
    input : variables parameter (optional, default = radiances)

    returns dictionary with data and meta
    """

    if not os.path.exists(filename):
        raise "Invalid Filepath"
    reader = SD(filename)
    aVariables = reader.datasets().keys()
    if variable==None:
        variable = 'radiances'
    elif variable in aVariables:
        pass
    else:
        raise "Invalid Variable Specified"

    data = reader.select(variable).get()
    #data = np.array(data)
    allz.append(data)
    outDict = {'Variable':variable,'filename':filename.split('/')[-1],'data':data}
    return outDict,np.vstack(allz)

Это определение get_lat_lon:

def get_lat_lon(path):
    allx = []
    ally = []
    reader = SD(path)
    lat = reader.select('Latitude').get()
    lon = reader.select('Longitude').get()    
    x,y = Proj(lon,lat)
    x /= 1000.0
    y /= 1000.0

    allx.append(x)
    ally.append(y)
    return np.vstack(allx),np.vstack(ally)

Это def interp_knn (который является приблизительным ближайшим соседом ANN)

def interp_knn(data, z, p):
    print "building kdtree"
    k = ann.kdtree(data)
    print "kdtree lookup..."
    ind, dist = k.knn(p, 1)
    print "done"
    img = z[ind[:,0]]
    img[dist[:,0] > 15] = N.NaN
    return img

и ошибка:

Traceback (most recent call last):
File "....\read_HDF5.py", line 166, in <module>
z.ravel(), np.column_stack((X.ravel(), Y.ravel())))
File "C:\Python27\lib\site-packages\numpy\lib\shape_base.py", line 296, in column_stack
return _nx.concatenate(arrays,1)
MemoryError

Значит, стек столбцов дает мне эту ошибку? Если это проблема, что я должен сделать, чтобы решить ее? Пожалуйста, дайте мне немного света.

---

EDIT:

Я набрал эти строки, чтобы распечатать каждое значение

print "x:",x
print "x.shape:",x.shape
print "y:",y
print "y.shape:",y.shape
print "X:",X
print "X.shape",X.shape
print "Y:",Y
print "Y.shape",Y.shape
print "x0:",x0    
print "xn:",xn    
print "y0:",y0    
print "yn:",yn

и я получил такой результат:

x: [[ 10424.20322635  10454.76060099  10485.45730949 ..., -12968.67726035
-12685.76602721 -12375.06502138]
[ 10382.59291927  10412.4034849   10442.35640928 ..., -12992.35321415
-12700.8632597  -12380.48805381]
[ 10340.74366218  10369.79366321  10398.98895233 ..., -13017.45507334
-12716.86098332 -12386.19350493]
..., 
[  5327.05493943   5275.15394042   5223.90854331 ...,   1918.57476975
1821.32106295   1717.34665908]
[  5303.06157859   5251.14693111   5199.89936454 ...,   1914.50352498
1818.19581363   1715.23546366]
[  5280.12577523   5226.55972784   5176.11746996 ...,   1910.4792526
1815.09866674   1714.77978295]]
x.shape: (135, 90)
y: [[ 8049.59989276  8099.28303285  8147.42741851 ...,  9925.58168202
9933.46845934  9937.89861612]
[ 8056.91586464  8106.78261584  8155.11136874 ...,  9953.01973235
9961.14109569  9965.68870206]
[ 8064.04624932  8114.09204498  8162.60060337 ...,  9980.50394667
9988.87543224  9993.54921283]
..., 
[ 7258.03197692  7292.42166577  7325.40914928 ...,  8225.26655004
8228.18675519  8230.16218915]
[ 7242.59306102  7276.75919255  7309.52794297 ...,  8201.49165135
8204.39528226  8206.36728948]
[ 7226.54007095  7261.56601577  7293.59601515 ...,  8177.75663252
8180.64399766  8182.58727191]]
y.shape: (135, 90)
X: [[-14149 -14148 -14147 ...,  14166  14167  14168]
[-14149 -14148 -14147 ...,  14166  14167  14168]
[-14149 -14148 -14147 ...,  14166  14167  14168]
..., 
[-14149 -14148 -14147 ...,  14166  14167  14168]
[-14149 -14148 -14147 ...,  14166  14167  14168]
[-14149 -14148 -14147 ...,  14166  14167  14168]]
X.shape (3635, 28318)
Y: [[ 7227  7227  7227 ...,  7227  7227  7227]
[ 7228  7228  7228 ...,  7228  7228  7228]
[ 7229  7229  7229 ...,  7229  7229  7229]
..., 
[10859 10859 10859 ..., 10859 10859 10859]
[10860 10860 10860 ..., 10860 10860 10860]
[10861 10861 10861 ..., 10861 10861 10861]]
Y.shape (3635, 28318)
x0: -14149
xn: 14168
y0: 7227
yn: 10861

Answer 1

Я согласен с yosukesabai в том, что вы должны распечатать x и y, но я думаю, что вы также усложняете ситуацию.Возможно, я не понимаю код, но кажется, что вы конвертируете все широтные координаты в км, а затем конвертируете латунные векторы из get_lat_lon в матрицу, а затем обратно в векторы.Я не думаю, что вам нужно это делать, по крайней мере, не для стандартной функции scipy kdtree.

Вот класс, который преобразует позиции в векторе lat-lon в соответствующие позиции i, j в сетке, где ячейки сетки имеют свои позиции, определенные матрицами lat-lon.Похоже, это то, что вы хотите.

Функция ll22ij вызывается с латентными векторами, соответствующими вашим данным.Затем вы можете использовать возвращенные пары ij для поиска значений в изображении с помощью img_matrix [ivec, jvec].

import numpy as np
import pylab as pl
from scipy.spatial import KDTree, cKDTree

import pycdf
from pyhdf.SD import SD,SDC

class SCB:
    def __init__(self, datadir="/projData/jplSCB/ROMS/",ijarea=[],
             lat1=None,lat2=None,lon1=None,lon2=None):
        self.i1 = 0     #
        self.i2 = 111   # Size of the grid.
        self.j1 = 0     # 
        self.j2 = 211   #
        self.datadir = datadir
        g = SD(datadir + '/scb_das_grid.nc', SDC.READ)
        self.lat = g.select('lat')[:]
        self.lon = g.select('lon')[:]-360
        self.llon,self.llat = np.meshgrid(self.lon,self.lat)

    def add_ij(self):
        i1=self.i1; i2=self.i2;j1=self.j1; j2=self.j2
        self.jmat,self.imat = np.meshgrid(np.arange(self.j2-self.j1),
                                          np.arange(self.i2-self.i1))
        self.ijvec = np.vstack((np.ravel(self.imat),np.ravel(self.jmat))).T

    def add_kd(self):
        self.kd = cKDTree(list(np.vstack((np.ravel(self.llon),
                                          np.ravel(self.llat))).T))
    def ll2ij(self,lon,lat,nei=1):
        if not hasattr(self,'kd'):
            self.add_kd()
            self.add_ij()
        dist,ij = self.kd.query(list(np.vstack((lon,lat)).T),nei)
        return self.ijvec[ij][:,0],self.ijvec[ij][:,1]

Причина для использования тычинок для add_kd и add_ij заключается в том, что генерацияkd-экземпляр для больших матриц.Я генерирую его только один раз, а затем повторно использую для разных наборов данных.основная концепция такова:

1. add_kd: cKDTree (или KDTree) инициируется с длинным списком длинных пар по лат (одна пара для каждой ячейки сетки).Эти пары генерируются путем выравнивания матриц lat и lon.

2. add_ij: две матрицы, состоящие из позиций i и j, сглаживаются так же, как матрицы lat и lon.

3. theвекторы со значениями lat и lon для наблюдений отправляются в функцию kd.query, и возвращается вектор с позициями ближайших пар.

Предположим, что следующая сетка, состоящая изиз трех матриц: координаты широты, долготы и данные:

---Lon---       ---Lat---    ---Data---   
12 13 14        30 30 30     5  8  3 
12 13 14        29 29 29     6  9  7
12 13 14        28 28 28     1  2  4

У нас есть наблюдения в следующей позиции по широте:

obs1: 12.2; 29.1
obs2: 13.4; 28.7

cKDtree будет запущено со следующей латпары lon:

12 28
13 28
14 28
12 29
13 29
14 29
12 30
13 30 
14 30 

и соответствующие пары ij будут

0  0
1  0
2  0
0  1
1  1
2  1
0  2
1  2
2  2

kd.query вернет

3 и 4,

, чтоэто позиции для координатных пар лат, которые наиболее близки к позициям наблюдений.Эти позиции также одинаковы в парах ij, что приводит к:

---Obs---         Grid
12.2; 29.1   ->   i=0, j=1
13.4; 28.7   ->   i=1, j=1

Поскольку сетка имеет следующие значения:

       5 8 3
vals = 6 9 7
       1 2 4

Теперь вы можете использовать vals [ivec,jvec], где ivec = [0,1] и jvec = [1,1], чтобы получить значения сетки, наиболее близкие к наблюдениям.ivec и jvec - это выход из ll2ij.

Answer 2

Мне кажется, что x0, xn, y0, yn - это проекционная координата в километрах. Затем вы построили X и Y, используя сетку с arange (x0, xn + 1), arange (y0, yn + 1). для каждого из этих диапазонов подразумевается разрешение в 1 км, поскольку размер шага диапазона равен 1, если не указано иное. Это то, что вы хотите? Это может быть огромный массив, например если я покрываю один континент с разрешением 1 км.

поэтому я предлагаю распечатку с проверкой реальности x, y и X, Y, посмотрите, выглядят ли они разумно.

EDIT

Посмотрев, что такое x, y и прочитав другие библиотеки Qs plus, которые вы использовали, я пришел к следующей версии. Я не могу проверить, потому что у меня нет данных modis. Так что дайте мне знать, если это не сработает, я собираюсь забрать то, что я пишу здесь, в этом случае. Попробуйте, пока ждете, пока решение brorfred сработает для вас.

if __name__=="__main__":

    filename = ... ( the hdf file I have)
    cumData, z = readAIRS_L1_VIS(filename)

    x, y = get_lat_lon(filename)

    # extent that satellite data covers
    x0, xn = x.min(), x.max()
    y0, yn = y.min(), y.max()

    # center point of data
    xo, yo = .5 * (x0+xn), .5*(y0+yn)

    # resolution of output grid, in km
    resolution = 20

    # ncol/nrow of image array
    ncol = int((xn - x0) / resolution) + 1
    nrow = int((yn - y0) / resolution) + 1

    # lower left corner of image array on projected coord
    p0 = xo - resolution * (ncol-1) * .5
    q0 = yo - resolution * (nrow-1) * .5

    # x,y coordinate of colomns and rows of image array on proj coord
    p = p0 + np.arange(ncol) * resolution
    q = q0 + np.arange(nrow) * resolution

    # x,y coordiate of all grid point of image array on projected coord
    X, Y = np.meshgrid(p, q)

    img = interp_knn(np.column_stack((x.ravel(), y.ravel())),
            z.ravel(), np.column_stack((X.ravel(), Y.ravel())))
    img.shape = (nrow, ncol)