GIS/RS knowledge base

集計など

  • ピボットテーブル的利用法
    pandas(ライブラリ)をインストールしておく
    Rやエクセルでやってもいいわけですが流した処理の集計結果を出力とかしたいのでpythonでやる
    (やってて思ったけどとても"R"っぽい)
 # ライブラリのインポート};~ 
import pandas as pd
 # CSV(出力結果)をデータフレームとして読み込み
df = pd.read_csv("C:\\temp\\data.csv")
 # とりあえず年と年毎に観測回数をカウントしてみる(NOは個体についてるID)
pd.crosstab([df2006.YEAR,df2006.MONTH], df2006.NO)
 # 月毎にグループ化(データフレーム全体)
grouped = df.groupby('MONTH')
 # 観測値(ここでは標高)の月毎平均を算出
grouped['elev'].mean()

データ操作関連

  • ファイル名の置換
    禁則文字を除去したいとき
    投影法などの情報をファイル名に付与したいとき
# os モジュールのインポート
import os
# サブディレクトリを含むファイル名を返す
pathiter = (os.path.join(root, filename)
     for root, _, filenames in os.walk("C:\\temp\\")
     for filename in filenames
)
# 特定の文字列を置換
for path in pathiter:
     # ハイフンをアンダーバーに変換指定
     newname =  path.replace('-', '_')
     if newname != path:
         os.rename(path,newname)
  • ファイルリストの取得
    • パスも含んだリスト(拡張子による指定も)
      import glob
      files_dir = glob.glob("C:\\temp\\*.csv")
      for file in files_dir:
          print file
    • 単純なファイルリスト
      import os
      files = os.listdir("C:\\temp\\")
      for file in files:
          print file
  • fnmatchをつかってもよい
    import os
    import fnmatch
    
    # データが保存されているパスをセット
    dir_p = "C:¥¥temp_ctr¥¥"
    # 拡張子を指定(*.shp)
    ext = "*.shp"
    # *.shpのファイル名を抽出しリストにする
    shplist = fnmatch.filter(os.listdir(dir_p), ext)
  • ex 上記のリストから特定の位置の文字列を取り出すとか
    for x in files_dir:
        start = len(x) - 12
        end = len(x) -4
        dftest = "df" + str(x[start:end])
        print dftest
  • ばらばらのDBFファイルをまとめてひとつのCSVへ変換
    dbfpy パッケージをインストールしておく
    # モジュールを呼び出す
    import csv
    from dbfpy import dbf
    import sys
    import glob
    
    # 出力CSVファイル
    fn_csv = "C:\\temp_c\\all_table.csv"
    # DBFファイルが格納されている場所
    path = "C:\\temp_c"
    
    # dbf ファイル名のリストを作成
    # 末尾に_zsと名前がついているdbfファイルのみを抽出 
    list_db = glob.glob(path + "\\*_zs.dbf")
    
    ## リストの最初のファイルをもとにCSVのひな形(列名など)を作る
    first_db = list_db[0:1]
    str_first_db = ",".join(map(str,first_db))
    init_db = dbf.Dbf(str_first_db)
    out_csv = csv.writer(open(fn_csv, 'wb'))
    
    # ひな形から列名リストを作りCSVファイルに書き込む
    names = []
    for field in init_db.header.fields:
        names.append(field.name)
    out_csv.writerow(names)
    
    # その他全てのdbfファイルからレコードのみをCSVファイルに書き込む
    for fn_db in list_db:
        in_db = dbf.Dbf(fn_db)
        for rec in in_db:
            out_csv.writerow(rec.fieldData)
        in_db.close()
        # 作業プロセスを表示
        print fn_db + "was converted"
    
    print "finsh all" 
  • 圧縮ファイル(tar.gz)の解凍
    # tarfileモジュールのインポート
    import tarfile
    # ファイルを指定する
    comp_f = tarfile.open("C:/temp/xxxx.tar.gz") 
    # 出力先を指定して解凍
    comp_f.extractall("C:/temp/")
  • 圧縮ファイル(zip)はこんな感じ
    import zipfile
    comp_f = zipfile.ZipFile("C:/temp/xxxx.zip")
    comp_f.extractall("C:/temp/")

pythonでベクタデータをちょす

pyshpをつかう

  • 属性によるデータのサブセット(属性によるフィーチャの選択→エクスポート)
# モジュールのインポート
import shapefile
import os
import fnmatch

# データが保存されているパスをセット
dir_p = "C:¥¥temp_ctr¥¥"
# 拡張子を指定(*.shp)
ext = "*.shp"
# *.shpのファイル名を抽出しリストにする
shplist = fnmatch.filter(os.listdir(dir_p), ext)

# データのサブセット処理
try:
    for i in shplist:        
        # データの読み込み
        r = shapefile.Reader(dir_p + str(i))
        # 書き込みインスタンスを定義(ポリライン)
        w = shapefile.Writer(shapeType=shapefile.POLYLINE)
        # 属性フィールドをコピー
        w.fields = list(r.fields)
        # 空のリストをつくる 
        sel_list = [] 
        for rec in enumerate(r.records()):
            # 100の倍数の標高値(8列目の属性フィールドに格納されている)を持つ等高線を抽出し書き込む
            if rec[1][7] % 100 == 0:
                sel_list.append(rec) 
        # ジオメトリと属性値を書き込む
        for rec in sel_list:
            w._shapes.append(r.shape(rec[0]))
            w.records.append(rec[1])
        # shapefileとして保存する(ファイル名の先頭にctr100mを付与)
        w.save(dir_p + "ctr100m" + str(i))
except:
    print "error!"
print "finish all"

Arcpyはじめのいっぽ

こんなときにはArcpyを使って処理するとべんりです。

  • バッチ処理をしたいとき
  • GUIで処理しようとしたら意味不明なエラーが出てきたとき
  • マシンのスペックが低くてGUIを起動するだけで精いっぱいなとき
    結構簡単なのでお試しください。
  • python用のエディタを用意するとよいです(Pythonインストール時に付いてくるIDLEでもよいです)
    テキストエディタでもよいですがpython用のエディタはインデント(pythonでは重要です)や関数などをわかりやすく整理・表示することができます。

pythonのおやくそく

  • 見やすいコードが好ましいです。リストなどが長くなる場合は適宜改行を入れるとよいです。
  • 括弧の中は改行できる
    	for x in [111, 222, 333, 444, 
    		555, 666, 777, 888, 999]:
  • 文字列の改行にはトリプルクォーテーションを使う
    """geometric
    correction"""
    みたいな感じ

各種処理のサンプル

シンプルな処理をしてみる

	# arcpyをインポートする
	import arcpy
	# 必要なエクステンションがあるかどうかをチェックする
	# spatial analyst = "spatial"
	# 3D analyst = "3D"
	arcpy.CheckOutExtension("spatial")
	# 入出力する変数を定義する
	# パスとファイル名はダブルコーテーションで囲う
	# パス区切りは\\で行う
	hexagons_shp = "D:\\test\\hexagons.shp"
	dem50gsi_j2ku54_img = "D:\\test\\dem50gsi_j2ku54.img"
	output_table = "D:\\temp\\zst_slope_02.dbf"
	# 処理コマンド中に上記変数を入力する
	# 上記の変数以外の設定値はダブルコーテーションで囲って指定する
	# どのような変数が必須かはHELPを参照する
	arcpy.gp.ZonalStatisticsAsTable_sa(hexagons_shp, "Unique_ID", dem50gsi_j2ku54_img, output_table, "NODATA", "ALL")
	# 処理が終了したらfinishという文字を表示させる
	print "finish"

たくさんのファイルをあつめてマージしてみる

  • ワイルドカード(*)を用いて複数ファイルを指定して実行する手順
	# globモジュールをインポートする
	import glob
	# arcpyをインポートする
	import arcpy
	# 出力ファイルのパスと名前を指定
	outputdata = "C:\\temp\\ps_merge_all.shp"
	# ワイルドカードを使用して入力データリストを指定
	files = glob.glob('C:\\temp\\*.shp')
	# Process: マージ (Merge)
	arcpy.Merge_management(files, outputdata)
	# 処理終了時にfinishと出力する
	print "finish"

バッチ処理をしてみる

ここでの例:国土基本メッシュ単位に分割されたレーザー測量データ(カンマ区切りテキスト)を読み込んでESRI shape file 形式に変換する

  • for文を使って複数ファイルに同じ処理を実施する
    	# arcpyをインポートする
    	import arcpy
    	# 入力および出力ファイルのディレクトリを指定
    	input_dir = "D:\\last_p"
    	output_dir = "D:\\last_p_out"
    	# 以下より処理内容
    	try:
    		for x in [111, 222, 333, 444, 555, 666, 777, 888, 999]:
            		#変数の宣言
            		fc = input_dir + "\\13lf" + str(x) + "_l_org_h.txt"
            		fc2 = str(x) + "_la"
            		fc3 = output_dir + "\\13lf" + str(x) + "_p_j2kxy13.shp"
            		# Process: XY イベント レイヤの作成(Make XY Event Layer
            		arcpy.MakeXYEventLayer_management(fc, "x", "y", fc2, "", "z")
            		# Process: フィーチャのコピー(Copy Features)
            		arcpy.CopyFeatures_management(fc2, fc3, "", "0", "0", "0")
            		# 処理が1サイクル正常に終了したら出力ファイル名を表示させる
            		print str(x) + "_p_j2kxy13.shp"
    	# どこかで処理が停止した場合 "error" と表示させる
    	except:
    	    print "error"  
    	# すべての処理が正常に終了したら"finish"と表示させる
    	print "finish"

入力データの生成など

250m間隔で最大60kmまでテキストを出力
(ex.250 500 750 1000・・・・)

# process for defining break value
# define minimum number
min_v = 250
# define maximum number
max_v = 60000
# define output text file
out_put = "C:\\temp\\output.txt"

# define range value
ini_v = min_v + min_v
end_v = max_v + min_v
int_v = min_v
# define break value
break_v = min_v

# set repeat number and interval
for x in range(ini_v,end_v,int_v):
    break_v = str(break_v) + " " + str(x)

# process text file output
f = open(out_put, 'w')
f.write(break_v)
f.close()
print break_v
print "finish"

フォルダ内すべてのshapefileにファイル名の入った属性フィールドを追加

import arcpy
import glob
import os.path
# set target directory
T_path = "C:\\temp\\"
# get file list as fullpath
listgps = glob.glob(T_path + "*.shp")
for f in listgps:
    # Process: Add Field
    arcpy.AddField_management(f, "orig_name", "TEXT", "", "", "100", "", "NULLABLE", "NON_REQUIRED", "")
    # get file name from fullpath
    f2 = os.path.basename(f)
    # set formula into a specific format
    f3 = "\"" + f2 + "\""
    # Process: Calculate Field
    arcpy.CalculateField_management(f, "orig_name", f3, "PYTHON", "")
    
print "finish"

for文に挿入するためCSVの列をリストに変換する

  • 特定のデータのみを選択して処理したい
  • 選択するファイルの名称やIDがスプレッドシートに入力されている
    というケースに使います
# 各モジュールをインポートする

import csv
import sys

# CSVファイルのフルパスを指定する

csvfile=open("C:\\temp\\list.csv")

# 各行のデータおよび各行の列数を格納するリストをそれぞれ宣言する

rowlist = []
colnum  = []

# 1行毎にリストにし列数をカウントする

for row in csv.reader(csvfile):
    rowlist.append(row)
    colnum.append(len(row))

# 行列変換 (各列の要素をリストにする)

collist = []
for colnum in range(max(colnum)):
    eachcollist = []
    for rownum in range(len(rowlist)):
        eachcollist.append(rowlist[rownum][colnum])
    collist.append(eachcollist)

# 行数(レコード数数)および列数の確認

print str(len(rowlist)) + " 行(レコード)" 
countcol = len(collist)

print countcol,"列",

# 指定した行をリストとして表示し確認
# ただし指定列数から1を引いた値を指定すること (1列目だったら0と指定)

print collist[0]

あとは for x in collist[0]: として繰り返しの処理を行う

polyline 各フィーチャの方位角を属性フィールドに追加する処理

→Arcgis10のツールボックス用です(そのように作るよう頼まれたため)
→各フィーチャは一直線という前提で書いています(フィーチャの始終点の座標から方位角を求めている)

import arcpy
import math

# Script arguments
input_table = arcpy.GetParameterAsText(0)

# local variables

infc = input_table

# input feature class

arcpy.AddField_management(input_table, "azimuth_a", "DOUBLE", "10", "5", "", "", "NULLABLE", "")

# Identify the geometry field
#
desc = arcpy.Describe(infc)
shapefieldname = desc.ShapeFieldName

# Create search cursor
#
rows = arcpy.UpdateCursor(infc)

# Enter for loop for each feature/row
#
for row in rows:
    # Create the geometry object
    #
    feat = row.getValue(shapefieldname)

    #Set start point
    startpt = feat.firstPoint

    #Set Start coordinates
    startx = startpt.X
    starty = startpt.Y

    #Set end point
    endpt = feat.lastPoint

    #Set End coordinates
    endx = endpt.X
    endy = endpt.Y

    radian = math.atan((endx - startx)/(endy - starty))
    degree = radian * 180 / math.pi
    row.azimuth_a = degree
    rows.updateRow(row)

# delete cursor and row objects to remove lock on the data
del row
del rows

GRASSのコマンドを pythonで動かす

vect = 'vector_data_sample'
    for x in ("TS199504","TS199505","TS199506","TS199507","TS199508","TS199509","TS199510","TS199511","TS199512","TS199601","TS199605"):
        grass.run_command("v.to.rast", input = vect, output = str(x) + "_cs", use = "attr", column = str(x))

python でジオコード

pygeocoderequestsパッケージが必要
→Google Geocoding APIを使っている

# 日本語でgeocodeしてみる
# -*- coding: utf-8 -*-

from pygeocoder import Geocoder
results = Geocoder.geocode(u"北海道富良野市")
print(results[0].coordinates)

#reverse geocodeしてみる
r_result = Geocoder.reverse_geocode(43.3421398, 142.3832254)
print(r_result[0])

トップ   編集 凍結解除 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2014-08-01 (金) 10:38:40 (3045d)