11/4 メシアではないので
救いを説くことは善にならない，なると信じることも自分自身できない．
希望を断つことが役割かもしれない，その誹りを受ける日がいつかやって来るのだろうか．
2431日と口にする人が，今でもそれを追い続けるわたしを導いてくれているというのは，とても幸福なことなんだろうと，思う．
川の名前と橋の名前
地味に位置と名前を一致させるのって手間がかかるけど，ここなら一覧表示．橋の名前もあって便利．
地図上にポリゴンを描いて，プロットしてある点がどのポリゴンに所属するかを確かめたい
1. 地図上にポリゴンを描く
まずフィーチャを新規作成する．ArcCatalog（ArcGISフォルダにあるはず）を起動，「新規作成」か右クリックでシェイプファイルを作成．今回はポリゴンを描きたいので，フィーチャタイプはポリゴンを選択．空間参照も地図に応じて決める．
そうしたらArcMapを開いて今作ったシェイプファイルをレイヤとして読み込む．
次に，「カスタマイズ」タブからエディターを出す．エディタ―タブから編集を開始．すると右端にフィーチャ作成ウィンドウが出て，その下に作図ツールが出るので，作成したいものを選択して操作，ポリゴンを描く．
このときのフィーチャ作成ウィンドウには同じフォルダにあるかつ現在のマップで読み込んでいるものがすべて出る．
エディタ―タブから保存して，終了を選択して，おしまい．
2. 描いたポリゴンの内部にある点を選び出す
「選択」タブから空間検索を行う．なお，このとき，シェイプファイルとして保存されていないと候補として出てこない（つまりcsv読み込んでxyデータを表示しただけの「…イベント」は候補として出ない）ので注意すること．
※シェイプファイルとしての保存は，コンテンツ（左側のウィンドウ）上で右クリック→データ→データのエクスポートから可能．
ターゲットレイヤーに点群，ソースレイヤーにポリゴンを入れて，ターゲットレイヤーフィーチャの空間選択方法を「ソースレイヤーフィーチャに含まれる」にすればOK．
3. 選択した点のデータをファイルにする
選択した点が属するレイヤを右クリック，表示してエクスポートで選択した点のみにするだけ．簡単．
各点までの距離をArcGISで出す
これまでは最近点しか出してこなかったけど，もちろん全部の点ペアについて距離を出すこともできる．
今回は各出発点Oから，各エリア重心への距離を求めてみた．
>>> import arcpy
>>> from arcpy import env
>>> env.workspace = "C:/data/pointdistance.gdb"
>>> inFeatures = "ODエリア分け完了IDつきGIS用_O"
>>> nearFeatures = "エリアデータ重心点"
>>> outTable = "O_重心点距離"
>>> searchRadius = "10000000000000 Meters"
>>> arcpy.PointDistance_analysis(inFeatures, nearFeatures, outTable, searchRadius)
最近接と同じ感じなのでさらりと．
このinとnearを同じにしたら，同セット内のペアとかもできそう．「入力フィーチャと近接フィーチャが同じレコードであるときはその結果がスキップされ、距離が 0 単位のフィーチャは報告されません」だそうです．
ちなみに
重心の出し方
これと同じで，各点の単純な平均にしちゃった．ただし外れ値を除いた．
ポリゴンの重心なら
このへんを見てやることになりそう．
最近点を求めたい場合の補足
この方法で，結合を使ってもできるみたい．あとGISの謎距離問題についても解説があって，
「JGD2000の設定された本州部分の東京の公示地価のシェープファイルを、平面直角9系に投影する。これは空間結合の距離計算時に、結合元の座標系の単位が使われるからだ。結合元は今回の場合、「各公示地価からの距離」をはかりたいので、公示地価のシェープファイルの単位が使われることになる。
ここで緯度経度のまま空間結合を実行すると、小数点以下のやたらと小さい距離が計算されて、悩み深くなることになるのでご注意を。」
だそうです．なるほど…
緯度経度から距離を計算する
ArcGISとかを使わない，Rのコードです．十進法表記で投入するだけ．
#Xが経度，yが緯度
D <- function(x1, y1, x2, y2){
#ラジアン表記に直す
x1 <- x1*pi/180
x2 <- x2*pi/180
y1 <- y1*pi/180
y2 <- y2*pi/180
P <- (y1 + y2)/2
dP <- abs(y1- y2)
dR <- abs(x1- x2)
m <- 6334834/sqrt*1^ 3)
n <- 6377397/sqrt(1 - 0.006674 * sin(P) * sin(P))
D <- sqrt*2
D
}
式は
これを参考にしました．これをデータ読み込んで自動で吐き出すようにしたいなら
#データに適用
ODdata <- read.csv("C:/workspace_R/new.csv")
distance <- matrix(0, nrow=length(ODdata[,1]), ncol=1)
for(i in 1:length(ODdata[,1])){
distance[i,1] <- D(ODdata[i,6], ODdata[i,7], ODdata[i,9], ODdata[i, 10])
}unlistData <- unlist(distance)
write.csv(unlistData, "C:/workspace_R//OD距離.csv")
とかで．
おまけ
これと，これを読む時間がほしいな～そろそろまた論文も読みたいな～
