10月 | 2014 | giscience.jp

ディレクトリに格納されている複数の航空レーザ測量のデータ(xxxx.txt)を
さっくりESRI shapefile形式(ポイント)に変換したいという相談があったので
pythonで書いてみました。

カンマで区切られたテキストの2列目がx 3列目がy 2列目がz
となっていたのでそれにあわせて書いてあります。

ogr2ogr使えばいいんだけど
筋トレ的に書いたので気にしない気にしない。

pyshpというライブラリが必要です。
コマンドプロンプトから
pip install pyshp
または
easy_install pyshp
でインストールできます。

pip も easy_install もインストールされていない場合は
このあたりからバイナリ(pyshp-1.2.0.win-amd64-py2.7.exe)をダウンロードして
インストールしてください。

その後以下をテキストエディタ等で適宜名前(ex. lp2shp.py)をつけて
Python27\Lib
直下に保存してください。

# 10/9追記：余計な処理をしてメモリを無駄に消費していたので修正

import shapefile
import csv
import glob
import os.path
def lp2shp(fd):
    # ファイルリストの作成
    flist = glob.glob(fd + "*.txt")
    # ファイルリスト毎に繰り返し処理
    for i in flist:
        # テキストファイルの読み込み
        cobj = csv.reader(open(i, "rb"))
        # ESRI shapefile(point)をつくる
        w = shapefile.Writer(shapefile.POINT)
        # 属性フィールド elev(標高値)を float, 出力フィールド幅:15, 精度:5, で作成する
        w.field("elev", "F", 15, 5)
        # 属性フィールド flagを integer 出力フィールド幅:5 で作成する
        w.field("flag", "N", 5)
        # テキストの行毎の情報をもとにポイント作成        
        for j in cobj:
            # 2列目をx座標 3列目をy座標としてポイント作成
            w.point(float(j[1]),float(j[2]))
            # 最初の属性フィールド(elev) に テキストファイル 4列目の情報を入力
            # 2番目の属性フィールド(flag) に テキストファイル 5列目の情報を入力
            w.record(float(j[3]),int(j[4]))
        # テキストファイルのフルパス情報からファイル名のみ抽出
        os.path.basename(i)
        # ESRI Shapefileの保存(テキストファイルのファイル名に"_p"を追加した名称を付与)
        w.save(str(i)[0:-4] + "_p.shp")
    print "finished"

if __name__ == '__main__':
    print "this is code block"

import shapefile

import csv

import glob

import os.path

def lp2shp(fd):

# ファイルリストの作成

flist = glob.glob(fd + "*.txt")

# ファイルリスト毎に繰り返し処理

for i in flist:

# テキストファイルの読み込み

cobj = csv.reader(open(i, "rb"))

# ESRI shapefile(point)をつくる

w = shapefile.Writer(shapefile.POINT)

# 属性フィールド elev(標高値)を float, 出力フィールド幅:15, 精度:5, で作成する

w.field("elev", "F", 15, 5)

# 属性フィールド flagを integer 出力フィールド幅:5 で作成する

w.field("flag", "N", 5)

# テキストの行毎の情報をもとにポイント作成

for j in cobj:

# 2列目をx座標 3列目をy座標としてポイント作成

w.point(float(j[1]),float(j[2]))

# 最初の属性フィールド(elev) にテキストファイル 4列目の情報を入力

# 2番目の属性フィールド(flag) にテキストファイル 5列目の情報を入力

w.record(float(j[3]),int(j[4]))

# テキストファイルのフルパス情報からファイル名のみ抽出

os.path.basename(i)

# ESRI Shapefileの保存(テキストファイルのファイル名に"_p"を追加した名称を付与)

w.save(str(i)[0:-4] + "_p.shp")

print "finished"

if __name__ == '__main__':

print "this is code block"

その後、python(command　line)　等で

import lp2shp as l
l.lp2shp("C:\\temp\\")

1 2	import lp2shp as l l.lp2shp("C:\\temp\\")

とテキストデータの入っているディレクトリを指定して実行してください。
同じディレクトリ内に xxxx_p.shp というファイル名でshapefileが出力されます。

前回の記事の続きです。

以下のような自動撮影装置の撮影記録 CSVファイルがあって

ID_A,地点名,ファイル名,日,時,日時,種名,ID_B,雌雄,行動
1,A,A1.avi,2014/01/01,10:00,2014/01/01 10:00,エゾシカ,3,メス,採餌中
2,A,A2.avi,2014/01/01,10:05,2014/01/01 10:05,エゾシカ,3,メス,採餌中
3,A,A3.avi,2014/01/01,11:00,2014/01/01 11:00,エゾタヌキ,3,オス,採餌中

ID_A,地点名,ファイル名,日,時,日時,種名,ID_B,雌雄,行動

1,A,A1.avi,2014/01/01,10:00,2014/01/01 10:00,エゾシカ,3,メス,採餌中

2,A,A2.avi,2014/01/01,10:05,2014/01/01 10:05,エゾシカ,3,メス,採餌中

3,A,A3.avi,2014/01/01,11:00,2014/01/01 11:00,エゾタヌキ,3,オス,採餌中

X番目に撮影された動画とX+1番目に撮影された動画の撮影時刻を比較し、
差が指定された時間以内(今回は分単位)
かつ
同種同性の個体が撮影されていれば
両者をグループ化する(末尾の列に同じID(gid)を追加する)というスクリプトをpythonで書いてみました。
データ処理後以下のような出力結果が得られます。

ID_A,地点名,ファイル名,日,時,日時,種名,ID_B,雌雄,行動,gid
1,A,A1.avi,2014/01/01,10:00,2014/01/01 10:00,エゾシカ,3,メス,採餌中,1
2,A,A2.avi,2014/01/01,10:05,2014/01/01 10:05,エゾシカ,3,メス,採餌中,1
3,A,A3.avi,2014/01/01,11:00,2014/01/01 11:00,エゾタヌキ,3,オス,採餌中,3

ID_A,地点名,ファイル名,日,時,日時,種名,ID_B,雌雄,行動,gid

1,A,A1.avi,2014/01/01,10:00,2014/01/01 10:00,エゾシカ,3,メス,採餌中,1

2,A,A2.avi,2014/01/01,10:05,2014/01/01 10:05,エゾシカ,3,メス,採餌中,1

3,A,A3.avi,2014/01/01,11:00,2014/01/01 11:00,エゾタヌキ,3,オス,採餌中,3

以下の通りメモとして残しておきます。

# -*- coding: utf8 -*-
##　標準モジュールのインポート 
import datetime
import csv
import dateutil.parser
### 入力例
## fpath  "C:\\temp\\sample.csv"
## intval  10
def gdata(fpath, intval):
    # csv ファイルの読み取り
    cobj = csv.reader(open(fpath, "rb"))
    # csv ファイルを2次元の配列データに変換
    arry = [ v for v in cobj]
    # 撮影間隔をdatetime形式に変換
    td = datetime.timedelta(0,0,0,0,intval,0,0)
    # 配列データの1行目にgid列を追加   
    ini = arry[0]
    ini.append("gid")
    # 繰り返し処理の初期値(0は列名行のため1からはじめる)
    i = 1
    x = 1
    # レコードの末端まで繰り返し処理
    while i < len(arry):
        # 配列のi行目にiというgidを挿入
        a = arry[i]
        a.append(str(i)) 
        # 次の行(レコード)がなければここで処理終了
        if x + 1 == len(arry):
            break
        # 次の行があれば処理を続行
        x2 = x + 1
        # x行目とx+1行目の6列目のdatetime文字列をdatetime形式にパース
        dto1 = dateutil.parser.parse(arry[x][5])
        dto2 = dateutil.parser.parse(arry[x2][5])
        # 同種同性で前のものと10分以内に撮影された個体であればiというグループ化ID(gid)を付与する
        while dto2 - dto1 < td and arry[x2][6]==arry[x][6] and arry[x2][8]==arry[x][8]:
            x = x + 1
            a = arry[x]
            a.append(str(i))
            # 次の行(レコード)がなければここで処理終了
            if x + 1 == len(arry):
                break
            # 次の行があれば処理を続行
            x2 = x + 1
            # x行目とx+1行目の6列目のdatetime文字列をdatetime形式にパース
            dto1 = dateutil.parser.parse(arry[x][5])
            dto2 = dateutil.parser.parse(arry[x2][5])            
        # 次に作成されるグループのリスト内での順番を定義
        i = x + 1
        x = x + 1
    # 同じフォルダに _g をつけたCSVファイルを新規作成 (同じ名前のファイルがあれば上書きする)
    f = open(fpath[0:-4] + "_g.csv", 'wb')
    # 配列データをCSVファイルに変換
    csvWriter = csv.writer(f)
    csvWriter.writerows(arry)    
    # CSVファイル書き込み終了
    f.close()
    print "finish" 

if __name__ == '__main__':
    print "this is code block"

# -*- coding: utf8 -*-

##　標準モジュールのインポート

import datetime

import csv

import dateutil.parser

### 入力例

## fpath "C:\\temp\\sample.csv"

## intval 10

def gdata(fpath, intval):

# csv ファイルの読み取り

cobj = csv.reader(open(fpath, "rb"))

# csv ファイルを2次元の配列データに変換

arry = [ v for v in cobj]

# 撮影間隔をdatetime形式に変換

td = datetime.timedelta(0,0,0,0,intval,0,0)

# 配列データの1行目にgid列を追加

ini = arry[0]

ini.append("gid")

# 繰り返し処理の初期値(0は列名行のため1からはじめる)

i = 1

x = 1

# レコードの末端まで繰り返し処理

while i < len(arry):

# 配列のi行目にiというgidを挿入

a = arry[i]

a.append(str(i))

# 次の行(レコード)がなければここで処理終了

if x + 1 == len(arry):

break

# 次の行があれば処理を続行

x2 = x + 1

# x行目とx+1行目の6列目のdatetime文字列をdatetime形式にパース

dto1 = dateutil.parser.parse(arry[x][5])

dto2 = dateutil.parser.parse(arry[x2][5])

# 同種同性で前のものと10分以内に撮影された個体であればiというグループ化ID(gid)を付与する

while dto2 - dto1 < td and arry[x2][6]==arry[x][6] and arry[x2][8]==arry[x][8]:

x = x + 1

a = arry[x]

a.append(str(i))

# 次の行(レコード)がなければここで処理終了

if x + 1 == len(arry):

break

# 次の行があれば処理を続行

x2 = x + 1

# x行目とx+1行目の6列目のdatetime文字列をdatetime形式にパース

dto1 = dateutil.parser.parse(arry[x][5])

dto2 = dateutil.parser.parse(arry[x2][5])

# 次に作成されるグループのリスト内での順番を定義

i = x + 1

x = x + 1

# 同じフォルダに _g をつけたCSVファイルを新規作成 (同じ名前のファイルがあれば上書きする)

f = open(fpath[0:-4] + "_g.csv", 'wb')

# 配列データをCSVファイルに変換

csvWriter = csv.writer(f)

csvWriter.writerows(arry)

# CSVファイル書き込み終了

f.close()

print "finish"

if __name__ == '__main__':

print "this is code block"

giscience.jp

月別アーカイブ: 2014年10月

航空レーザ測量データをひたすら ESRI shapefile に変換するだけ

カメラトラップの撮影記録を撮影間隔&種でグループ化