Pythonで大きなデータセットから2次元の緯度、経度、データをサブセット化する（IMS snowcover）

Question

私は「シンプルな」目的を念頭に置き、達成することは非常に困難です。私は3つの2次元配列lat,lonおよびdataをすべて次元（24576、24576）で持っています。最初の2つはマップ上にプロットしようとしている変数dataの緯度と経度座標です。私はバイナリファイルとテキストファイルの組み合わせからこれらのデータをすべて読み込んでいます。そのため、前処理操作は実質的に不可能であり、Pythonスクリプトで行う必要があります。

配列の次元が与えられると、メモリの制約のため、地球上の小さな領域でベースマップ投影を選択してもデータを直接描画することは事実上不可能です。私はすでにそれを試して、basemap.contourfでプロットしようとするとメモリエラーが発生しました。

したがって、輪郭関数に渡す前に配列をサブセット化する必要があります。私は多くのことを試しましたが、何も動作していないようです。私の考えはlatsとlonsは、2-Dの配列座標です。この

lat_bnds, lon_bnds = [35, 50], [5, 20] 
condition=((lats > lat_bnds[0]) & (lats < lat_bnds[1])) & (lons > lon_bnds[0]) & (lons < lon_bnds[1]) 

ような何かをしました。これは、元のデータをマスクするために使用することができるが、サブセットをマスクすることができないlats（または同等にはlons）と同じ次元を持つ2次元ブール値配列を生成する。

で同じ条件を使用すると、使用できないアレイ形状の(2, 2564856)が生成されます。ここで問題となるのは、2-D配列のあらゆる点で満足される必要がある複数の条件があり、連続する矩形部分行列につながるという保証はないということです。しかし、それらのデータをプロットするために、私は実際にそれらを部分集合するか、またはそれらを詳しく説明する別の方法を見つける必要があります。

明らかなものがありませんか？エラーが発生することなく元のデータをプロットする他のスマートな方法はありますか？あなたが必要なファイルをダウンロードしたと仮定して、データを読み込むためのhttp://nsidc.org/data/G02156

簡単なスクリプト：データの

ソース

import numpy as np 
import pandas as pd 
from mpl_toolkits.basemap import Basemap # Import the Basemap toolkit 

with open('IMS1kmLats.24576x24576x1.double', 'rb') as f: 
data = np.fromfile(f, dtype='<d', count=24576*24576) 
lats = np.reshape(data, [24576, 24576], order='F') 

with open('IMS1kmLons.24576x24576x1.double', 'rb') as f: 
data = np.fromfile(f, dtype='<d', count=24576*24576) 
lons = np.reshape(data, [24576, 24576], order='F') 

widths=np.full((24576), 1, dtype=int).tolist() 
data=np.array(pd.read_fwf('ims2017312_1km_v1.3.asc', skiprows=30, 
widths=widths, lineterminator='\n', header=None)) 

Answer 1

私はあなたの配列を得た方法を確認していません元の配列のサブセットを必要とし、条件を満たす場合は、元のリストからスライスすることができます。

a = np.random.randint(10,50, size=(50,2)) # ---- generate coordinates 

array([[44, 11], 
     [40, 36], 
     [19, 26], 
     ..., 
     [33, 26], 
     [42, 12], 
     [15, 25]]) 

lats = a[:, 1] # ---- latitude is Y-axis 
lons = a[:, 0] # ---- longitude is X-axis 
lat_bnds, lon_bnds = [35, 50], [5, 20] # ---- your desired bounds 

condition =((lats > lat_bnds[0]) & (lats < lat_bnds[1])) & 
      (lons > lon_bnds[0]) & (lons < lon_bnds[1]) 

condition 
array([False, False, False, ..., False, False, False], dtype=bool) 

condition.shape # => (50,) 

a[condition]  # ---- slice the original array 

array([[11, 45], 
     [15, 40], 
     [15, 43], 
     [15, 49]]) 

あなたが行きたい方向にあなたを得るための希望。

Answer 2

PythonでIMS 1 kmの解像度データをプロットしようとしている間に誰かがこのポストを見つけた場合、うまくいくはずの貧弱な解決策があります。

プロットルーチンによってメモリエラーが発生しますが、配列はまだPythonに格納できます。だから、代わりにどこ機能を使用してサブセット化、私はちょうどのような

lat_subset=lat[imin:imax, jmin:jmax] 

を明示的にインデックスを使用しようとしましたし、その後持っているために地図投影を等高線図をやったり使用せずにplot.imshow（）で結果をプロットデータがどのように見えていたかのアイデア。これにより、関心のある領域が内部にあるようにインデックスのスパンを選択することができました。今私は、メモリエラーを起こさずに等高線図を作ることができました。

これらのデータをプロットする方法を示すノートブックがある小さなリポジトリを持っています：https://github.com/guidocioni/snow_ims/blob/mistral/plot_ims.ipynb。 座標ファイルはサイズを指定してダウンロードする必要がありますが、ファイルをオンラインで直接読むことができるという利点があります。