python - 是否有任何与在 numpy PYTHON 中使用“列堆栈”相同的方法

Question

我正在使用 Python 2.7。我的系统正在运行 Window Vista，32 位。

我有一段代码可以读取辐射度、纬度和经度以及一个图像文件（以 hdf 扩展名）。然后尝试执行近似最近邻并对其进行映射。但是当它试图做近似最近的邻居时，它给了我内存错误。

仅 hdf 文件就有 4.70 MB，看起来大小不算太大。

这是我的代码：

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)

这是 def 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

那么列堆栈是否给了我这个错误？如果这是问题所在，我应该怎么做才能解决它？请给我一些光。

编辑：

我输入了这些行以打印出每个值

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

score 2 · Accepted Answer

我同意 yosukesabai 的观点，即你应该打印出 x 和 y，但我认为你也让事情变得比必须做的更难。我可能不理解代码，但似乎您正在将所有经纬度坐标转换为 km，然后将经纬度向量从 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 的 if 雄蕊的原因是为大型矩阵生成 kd 实例的成本很高。我只生成一次，然后将其重用于不同的数据集。基本概念如下：

add_kd：cKDTree（或 KDTree）由一长串经纬对（每个网格单元一对）启动。这些对是通过展平 lat 和 lon 矩阵生成的。
add_ij：由 i 和 j 位置组成的两个矩阵以与 lat 和 lon 矩阵相同的方式展平。
具有观测值的 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 将使用以下 lat-lon 对启动：

相应的 ij 对将是

kd.query 将返回

3和4，

这是最接近观测位置的网格 lat-lon 对的位置。这些位置在 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

You can now use vals[ivec,jvec] where ivec=[0,1] and jvec=[1,1] to get the grid values that closest corresponds to the observations. ivec and jvec is the output from ll2ij.

score 1 · Accepted Answer

在我看来，x0、xn、y0、yn 是以公里为单位的投影坐标。然后，您使用带有 arange(x0,xn+1)、arange(y0,yn+1) 的网格网格构造 X 和 Y。这些范围中的每一个都隐含假设为 1km 分辨率，因为范围的步长为 1，除非您另有说明。这是你想要的吗？它可能是一个巨大的阵列，例如，如果我以 1 公里的分辨率覆盖一个大陆。

所以我建议现实检查打印 x, y 和 X, Y，看看它们是否合理。

编辑

在查看了 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)

python - 是否有任何与在 numpy PYTHON 中使用“列堆栈”相同的方法

2 回答 2

Related

Reference