![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
集計など †
# ライブラリのインポート};~
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)
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を使って処理するとべんりです。
pythonのおやくそく †
各種処理のサンプル †シンプルな処理をしてみる † # 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 形式に変換する
入力データの生成など †250m間隔で最大60kmまでテキストを出力 # 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の列をリストに変換する †
# 各モジュールをインポートする
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 でジオコード †→pygeocode と requestsパッケージが必要 # 日本語で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]) |
