有没有与numpy中“列堆叠”相同的函数？
我在用Python 2.7，系统是Windows 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)
这是定义的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)
这是定义的interp_knn函数（就是近似最近邻的算法）
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
2 个回答
1
我觉得x0、xn、y0、yn是以公里为单位的投影坐标。然后你用meshgrid和arange(x0,xn+1)、arange(y0,yn+1)构建了X和Y。每个arange的隐含假设是分辨率为1公里，因为arange的步长默认是1，除非你另行指定。这样做是你想要的吗？如果我以1公里的分辨率覆盖一个大陆，那可能会生成一个非常大的数组。
所以我建议你先打印一下x、y和X、Y，看看它们是否合理。
编辑
在了解了x和y的含义，以及你其他问题和使用的库后，我想出了以下版本。我无法测试，因为我没有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)
2
我同意yosukesabai的看法，你应该打印出x和y的值，但我觉得你把事情搞得比实际复杂了一些。我可能不太懂代码，但看起来你是把所有的经纬度坐标转换成公里，然后把从get_lat_lon得到的经纬度向量转换成矩阵，再转换回向量。我觉得你不需要这样做，至少对于标准的scipy kdtree函数来说，不需要。
这里有一个类，可以把经纬度向量中的位置转换成网格中的对应i,j位置，网格单元的位置是由经纬度矩阵定义的。这似乎正是你想要的。
函数ll22ij会用与你的数据对应的经纬度向量来调用。然后你可以用返回的i-j对在图像中查找值，比如用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）是用一长串经纬度对初始化的（每个网格单元一个对）。这些对是通过将经纬度矩阵展平生成的。
add_ij: 两个包含i和j位置的矩阵也以与经纬度矩阵相同的方式展平。
观察值的经纬度向量会被发送到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会用以下经纬度对初始化：
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，
这就是与观察值位置最近的网格经纬度对的位置。这些位置在i-j对中也是相同的，得出：
---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的输出。