ShapelyとBasemapによるGISデータ操作1

インストール

CentOS 7上のPython 3またはPython 2で地理情報データを扱うため，Shapelyとbasemapをcondaとpipでインストールします．Anacondaをインストール済であることを前提とします．
[cc]
> conda install gdal
> conda install fiona
> conda install basemap
> conda install shapely
> pip install descartes
[/cc]
basemapでgeos3.3.3-0が必要のようにエラーが出る場合は再度basemapをインストールしてgeosのバージョンを下げる必要があるようです．
[cc]
> conda install basemap
[/cc]
以降の解析のためｍ，以下のパッケージが読み込まれていることを前提とします．
[cc]
import fiona
import numpy as np
from pandas import Series, DataFrame
import pandas as pd
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
[/cc]

shapefileを読み込みDataFrameを作成

ArcGISのshapefileを扱うには，まずFionaで全体を読み込みます．その際，*.shx等の関連ファイルが同じディレクトリに存在している必要があります．
[cc]
In [1]: fname = 'TokyoBay5m_original'
fshp = fname + '.shp'
shp = fiona.open(fshp) # shapefileをFionaで読み込む
points = [c for c in shp] # 1行を要素とするリストを作る
[/cc]
pointsの要素がどのような構造をしているか，最初の要素を表示して調べます．
[cc]
In [2]: points[0]
Out[2]:
{'geometry': {'coordinates': (-19000.0, -90200.0), 'type': 'Point'},
'id': '0',
'properties': OrderedDict([('I', 6),
('J', 1),
('X', -19000),
('Y', -90200),
('DEPTH', 0.0),
('MASK', 0)]),
'type': 'Feature'}
[/cc]
pointsの要素は辞書であることが分かります．例えば，辞書のキー「properties」にある，I，J，X，Y，DEPTH，MASKの値を抽出してリストを作ります．さらに，この例では要素はすべて数値なので，高速化のため，NumPyのndarrayに変換しておきます．なお，forループを使いましたが，リスト内包表記でも作れます．ただし，pointsのサイズが大きいと，リスト内包では（おそらく）メモリ不足からうまく動かないことがありましたので，forループの方が安全でしょう．
[cc]
Mtmp = []
for n in range(len(points)):
Mtmp.append(list(points[n]['properties'].values()))
M = np.array(Mtmp) # Matrix Mtmpの要素をndarrayに変換してMatrix Mを作る
[/cc]
Matrix Mの中身を見てみましょう．
[cc]
In [3]: M
Out[3]:
array([[ 6.00000000e+00, 1.00000000e+00, -1.90000000e+04,
-9.02000000e+04, 0.00000000e+00, 0.00000000e+00],
[ 7.00000000e+00, 1.00000000e+00, -1.89500000e+04,
-9.02000000e+04, 0.00000000e+00, 0.00000000e+00],
[ 8.00000000e+00, 1.00000000e+00, -1.89000000e+04,
-9.02000000e+04, 0.00000000e+00, 0.00000000e+00],
...,
[ 9.33000000e+02, 7.80000000e+02, 2.76000000e+04,
-2.60500000e+04, 0.00000000e+00, 0.00000000e+00],
[ 9.34000000e+02, 7.80000000e+02, 2.76500000e+04,
-2.60500000e+04, 0.00000000e+00, 0.00000000e+00],
[ 9.35000000e+02, 7.80000000e+02, 2.77000000e+04,
-2.60500000e+04, 0.00000000e+00, 0.00000000e+00]])
[/cc]
得られたリストにキーを付けて（以下の例ではIをi，Jをjのように変えました）辞書のリストを再構築し，最後にpandasのDataFrameに変換します．リスト内包表記を使っています．
[cc]
i = [row[0] for row in M]; j = [row[1] for row in M]
x = [row[2] for row in M]; y = [row[3] for row in M]
depth = [row[4] for row in M]; mask = [row[5] for row in M]
data = {'i':i, 'j':j, 'x':x, 'y':y, 'Depth':depth, 'Mask':mask}
df = DataFrame(data)
[/cc]
最後にdfの中身を見てみましょう．実際はもっと長く表示されますが，省略しています．
後はpandasを使って効率的にデータのハンドリングができます．
[cc]
In [4]: df
Out[4]:
Depth i j Mask x y
0 0 6 1 0 -19000 -90200
1 0 7 1 0 -18950 -90200
2 0 8 1 0 -18900 -90200
3 0 9 1 0 -18850 -90200
... ... ... ... ... ... ...
1200537 0 933 780 0 27600 -26050
1200538 0 934 780 0 27650 -26050
1200539 0 935 780 0 27700 -26050

[1200540 rows x 6 columns]
[/cc]