0

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)

これは定義 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 です (近似された最近傍 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
4

2 に答える 2

2

x と y を出力する必要があるという yosukesabai に同意しますが、必要以上に物事を難しくしていると思います。コードが理解できないかもしれませんが、すべての緯度経度座標を km に変換し、緯度経度ベクトルを get_lat_lon から行列に変換してからベクトルに戻しているようです。少なくとも標準の scipy kdtree 関数では、そうする必要はないと思います。

以下は、緯度経度ベクトルの位置をグリッド内の対応する i,j 位置に変換するクラスで、グリッド セルの位置は緯度経度行列によって定義されます。これはあなたが望むもののようです。

関数 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 インスタンスを生成するのにコストがかかるためです。一度だけ生成してから、別のデータセットに再利用します。基本的なコンセプトは次のとおりです。

  1. add_kd: cKDTree (または KDTree) は、緯度と経度のペア (グリッド セルごとに 1 つのペア) の長いリストで開始されます。これらのペアは、lat 行列と lon 行列を平坦化することによって生成されます。

  2. add_ij: i と j の位置で構成される 2 つの行列は、緯度と経度の行列と同じ方法で平坦化されます。

  3. 観測値の緯度と経度の値を持つベクトルが kd.query 関数に送信され、最も近いペアの位置を持つベクトルが返されます。

緯度、経度の位置、およびデータの 3 つのマトリックスで構成される次のグリッドを想定してみましょう。

---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、

これは、観測の位置に最も近いグリッドの緯度と経度のペアの位置です。これらの位置は 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] where ivec=[0,1] および jvec=[1,1] を使用して、観測に最も近いグリッド値を取得できるようになりました。ivec と jvec は ll2ij からの出力です。

于 2011-10-30T02:01:59.703 に答える
1

x0、xn、y0、yn がキロメートル単位で投影されているように見えます。次に、arange(x0,xn+1)、arange(y0,yn+1) のメッシュグリッドを使用して X と Y を構築しました。特に指定しない限り、arange のステップ サイズは 1 であるため、これらの各 arange には 1 km の解像度の暗黙の仮定があります。これは、あなたの望むことですか？たとえば、1 km の解像度で 1 つの大陸をカバーする場合など、巨大な配列になる可能性があります。

したがって、x、y および X、Y を印刷して現実を確認することをお勧めします。それらが合理的に見えるかどうかを確認してください。

編集

x、yが何であるかを調べ、使用した他のQとライブラリを読んだ後、次のバージョンを思いつきました。modisデータがないため、テストできません。これが機能しない場合はお知らせください。その場合、ここに書いたことを撤回します。brofred のソリューションが機能するのを待つ間、これを試してください。

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)
于 2011-10-30T01:16:09.320 に答える