# ArcGISでの座標変換
## 概要
平面直角座標での津波計算を行うために，断層パラメータの緯度・経度を`ArcGIS`を用いて座標変換しました．
`ArcMap`を使ってマウスでいちいちやるのが面倒だったので,`arcpy`でスクリプト書いてやりました．
今回やったのは，断層の左下座標の変換です．
## 使用データ
- [Fujii,Satake モデル 断層パラメータ Ver8.0 ](http://iisee.kenken.go.jp/staff/fujii/OffTohokuPacific2011/tsunami_ja.html)
※事前にcsvファイルにしておく
## 環境
- ArcMap 10.5
## 手順
1.csvファイルをポイント( lyr )に変換
2.投影座標に変換( lyr -> shp )
3.XYデータをテーブルに追加
4.テーブル( shp )をcsvファイルに変換
## csvファイルをポイント( lyr )に変換
<details>
最初のSTEPだが，ここが一番苦労した...
まず，以下のサイトを参考にコードを次のように書いた．
[XYテーブル→ポイント(XY Table To Point)-ヘルプ｜ArcGIS](http://pro.arcgis.com/ja/pro-app/tool-reference/data-management/xy-table-to-point.htm)
```python
import arcpy
#
# setting
in_table = 'input.csv'
out_feature_class = 'out_table'
x_field = 'long'
y_field = 'lat'
coor_sys = arcpy.SpatialReference('WGS 1984')
#
# excute
arcpy.management.XYTableToPoint(in_table,out_feature_class, \
x_feild, y_field,coor_system)
```
各変数の説明は面倒なので，サイトを参考にしてください．
「案外簡単じゃん」と回してみたら，次のエラーが出た．
```
AttributeError: 'module' object has no attribute 'XYTableToPoint'
```
`XYTableToPoint`モジュールがない，だと...
いやいや，公式サイトに書いてあるんだからあるでしょうよ！？
しかし，ArcMapのPythonウインドウで確認したらなかった...
なんでないのよ，と思いながら最初からやり直し．
時間かかった結果，以下のサイトを発見.
[XYイベントレイヤーの作成 (Make XY Event Layer)-ヘルプ｜ArcGIS](http://desktop.arcgis.com/ja/arcmap/10.3/tools/data-management-toolbox/make-xy-event-layer.htm)
これを参考に以下のコード作成．
```python
import arcpy
#
# Setting
in_table = 'input.csv'
out_layer = 'out_table'
x_field = 'long'
y_field = 'lat'
coor_sys =arcpy.SpatialReference('WGS 1984')
save_layer = 'out_point.lyr'
#
# make layer
arcpy.MakeXYEventLayer_management(in_table, x_field, y_field, \
out_layer, coor_sys)
#
# save layer
arcpy.SaveToLayerFile_management(out_layer,save_layer)
```
今度は成功したコードなので，軽く説明を．
(x,y)=(long,lat)というポイントのレイヤーを作成したい．
`make layer`でそいつを作ります．
しかしこれは一時的なレイヤーで，保存しない限りどこのフォルダにも存在しないので，`save layer`で保存します．
`coor_sys`は座標系の定義をするもので，引数には上記コードのように文字を正しく与えるか，EPSGで定義されている数字を与えます．どっちを選ぶかは好みじゃないかと...
EPSGの説明は次を読んでください．[EPSGコードとは？｜空間情報クラブ](http://club.informatix.co.jp/?p=1225)
</details>
## 投影座標に変換( lyr -> shp )
<details>
このSTEP以降は無駄口少なめにいきます．
ここで行うことが最も重要ですね，簡単だったけど．
参考にしたのは次のサイト．[投影変換 (Project)-ヘルプ|ArcGIS](http://desktop.arcgis.com/ja/arcmap/10.3/tools/data-management-toolbox/project.htm)
作成したコードは以下の通り．
```python
import arcpy
#
in_dataset = 'out_point.lyr'
out_dataset = 'out_poject.shp'
out_coor_sys= arcpy.SpatialReference(6678) # JGD_2011_Japan_Zone_10
#
arcpy.Project_management(in_dataset,out_dataset,out_coor_sys)
```
inputデータには先ほど作成した`lyr`ファイルを指定．
outputデータには`shp`ファイルを指定．
ファイル名には変換した座標系の名前を入れるとわかりやすいかと．
また`shp`ファイルは色々とファイルを作成するので，フォルダに1つずつ作成しましょう．
</details>
## XYデータをテーブルに追加
<details>
続いて，変換したXY座標のデータを`shp`ファイルのテーブルに追加します．
参考にしたのは次のサイト．[XY 座標の追加 (Add XY Coordinates)-ヘルプ｜ArcGIS](https://pro.arcgis.com/ja/pro-app/tool-reference/data-management/add-xy-coordinates.htm)
```python
import arcpy
#
in_feature = 'out_project.shp'
out_feature = 'out_project_addxy.shp'
#
arcpy.Copy_management(in_feature,out_feature)
arcpy.AddXY_management(out_feature)
```
作った`shp`ファイルをコピーして，新しい`shp`ファイルのテーブルにXY座標を追加します．
コピーしなくてもいいですが，なんかあった時のために最初はコピーをおすすめします．
特に問題なければ，コピーせずにやってもいいです．`shp`ファイルは場所とるので．
その時はこの作業は，`import`含めて3行で終わります．引数も少ないし，簡単ですね．
</details>
## テーブル( shp )をcsvファイルに変換
<details>
最後のSTEPです．`csv`にして出力しましょう．
参考サイトは以下の2つ．
[テーブル ビューの作成 (Make Table View)-ヘルプ|ArcGIS](https://pro.arcgis.com/ja/pro-app/tool-reference/data-management/make-table-view.htm)
[テーブル→テーブル(Table to Table)-ヘルプ|ArcGIS](http://desktop.arcgis.com/ja/arcmap/10.3/tools/conversion-toolbox/table-to-table.htm)
```python
import arcpy
#
in_shp = 'out_project_addxy.shp'
temp_table = 'temp_table'
out_path = './'
out_name = 'output.csv'
#
arcpy.MakeTableView_management(in_shp, temp_table)
#
arcpy.TableToTable_conversion(temp_table, out_path, out_name)
```
以上のコードでは，XY座標を追加した`shp`ファイルから保存されないテーブル`temp_table`を作成し，さらにそれを`csv`に変換して出力している．
`shp`ファイルから直接は`csv`にできないようなので，このようにした．
出力された`csv`には`POINT_X`,`POINT_Y`の列があるので，それが変換された座標です．
</details>
## まとめ
<details>
今回はArcGISの`arcpy`を用いたポイントデータの座標変換を行った．
`arcpy`を使うのは初めてだったので，調べるのに苦労したけれど，コードは意外と簡単にできましたね．
マウスを使わなくてもいいのは，すごく快適でした．
STEP毎に分けて行ったので，長い説明になったけど一括のコードにすると，もっと早く分かりやすいコードになります．
それも書いたので，最後に添えます．
説明は特につけないので，わかんないとこは調べてください．
ではでは．
```python
import arcpy
import os
#
# ----- Setting -----
# name of input and output csv
in_csv = 'input.csv'
out_csv = 'output.csv'
#
# coordinate system
in_coor_sys = arcpy.SpatialReference('WGS 1984')
out_coor_sys = arcpy.SpatialReference('Tokyo_UTM_143.prj')
#
# ----- excute (please no change) -----
# make temporary directory
os.mkdir('temp')
temp_d = 'temp/'
# csv to point layer
x_field = 'long'
y_field = 'lat'
temp_layer = 'point_layer'
#
print('now processing : csv to point layer')
arcpy.MakeXYEventLayer_management(in_csv,x_field, y_field,\
temp_layer, in_coor_sys)
#
# projection in_coor_sys to out_coor_sys
temp_shape = temp_d + 'point.shp'
#
print('now processing : projection from ', in_coor_sys, \
' to ', out_coor_sys )
arcpy.Project_management(temp_layer, temp_shape, out_coor_sys)
#
# add XY data to shape file
print('now processing : add XY data to shape file ')
arcpy.AddXY_management(temp_shape)
#
# point shape to out csv
temp_table = 'temp_table'
out_path = './'
#
print('now processing : point shape to output csv')
arcpy.MakeTableView_management(temp_shape,temp_table)
arcpy.TableToTable_conversion(temp_table, out_path, out_csv)
#
# remove temporary directory
#os.rmdir('temp')
```
</details>
###### tags : `ArcGIS`
Sign in with Wallet
Connect another wallet