経度緯度をもつGeoDataFrameを平面直角座標系（JGD2011）に変換するPythonツールtoflatを作りました

日本の経度緯度の位置情報を扱っていると、平面直角座標系に変換したい場面がよくあります。たとえば２点間の距離を知りたいような場合に、変換作業を行います。

しかし、次の地図に示すように、平面直角座標系は19個あるため、各位置の番号を覚えるのは（私には）難しいことでした。

![スクリーンショット 2023-11-10 093731.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/85433/a327f9d9-fbbb-62ae-4cd5-1f8d73591839.png)


というわけで、この作業のたびにウェブでどの番号かを調べて割り当て作業をしていたのですが、これって非効率だなぁと思うようになりました。

そこで、GeoDataFrameに持ったデータに対して、最も近い平面直角座標系を割り当ててくれるPythonのツールtoflatを作成して、PYPIに公開してみました。

https://github.com/mazarimono/toflat

## 使い方
### インストール
まずはtoflatをインストールしてください。

```
pip install toflat
```

### toflatの動作概要

GeoDataFrameの位置情報からcentroidを作り、それに近い各座標の原点を探す。

### 使ってみる

次に京都の名所の位置情報のGeoDataFrameを作成します。

```python
from toflat import toflat
import geopandas as gpd
from shapely.geometry import Point

point1 = Point(135.758, 34.985)
point2 = Point(135.729, 35.039)
point3 = Point(135.746, 35.066)
point4 = Point(135.798, 35.027)

data = {
    'name': ['Kyoto Station', 'Kinkaku temple', 'MK bowl', 'Ginaku_temple'],
    'geometry': [point1, point2, point3, point4]
}
gdf = gpd.GeoDataFrame(data, geometry='geometry').set_crs('EPSG:4326')
```

上で作成したGeoDataFrame: gdfは次の画像のようになります。

![スクリーンショット 2023-11-10 094322.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/85433/937877c4-40ec-e356-b3e1-4ede0a163b0b.png)


CRSも確認しておきましょう。
```python
gdf.crs
```
```
<Geographic 2D CRS: EPSG:4326>
Name: WGS 84
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: World.
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984 ensemble
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich
```

はい。次にtoflatを使って、平面直角座標系に変換します。使い方は簡単で、ToFlatインスタンスにCRSを変換したいデータを渡して、to_flatすると、平面直角座標系に変換したデータが作られます。

```python
tf = toflat.ToFlat(gdf)
tf_gdf = tf.to_flat()
print(tf_gdf)
print(tf_gdf.crs)
```

実行結果は次の通りです。

![スクリーンショット 2023-11-10 095020.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/85433/73f8e262-e157-d1af-e67b-4446bcc04234.png)


### 直線距離を求めてみる

京都駅からの距離を経度緯度のままのデータ（gdf）と、toflatで作ったデータ(tf_gdf)を使って求めてみます。経度緯度だとよくわからない情報ですが、平面直角座標系に変換したらメートルで距離が出せます。

```python
kyoto_station = gdf.loc[0, 'geometry']
gdf['distance'] = gdf['geometry'].map(lambda x: kyoto_station.distance(x))
gdf
```
![スクリーンショット 2023-11-10 095608.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/85433/e8c79273-f265-0ce9-9732-3387942ff064.png)


```python
kyoto_station = tf_gdf.loc[0, 'geometry']
tf_gdf['distance'] = tf_gdf['geometry'].map(lambda x: kyoto_station.distance(x))
tf_gdf
```
![スクリーンショット 2023-11-10 095643.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/85433/4888591a-5d94-8f0a-8555-47039ea89397.png)

## まとめ
というわけで、平面直角座標系に変換するツールを作成して公開しました。まだいろいろな場面で使ってないので、いろんなケースでどうなるか分かってないので、その辺りのテストを今後は書きたいと思っています。

CRSの戻し方ですが、つぎのようにGeoPandasにあるto_crsメソッドを使ってできます。

```
tf_gdf4326 = tf_gdf.to_crs('EPSG:4326')
```

toflatを試してもらって、変な動作が見つかったら、記事のコメントissuesにいただけますと幸いです（スターも大歓迎です）。

https://github.com/mazarimono/toflat

記事がちょっとでも参考になったら「いいね、筆者のフォロー」をしていただけますと幸いです。

**お読みいただきありがとうございました。**

今回のコード（Google Colab）

https://colab.research.google.com/github/mazarimono/toflat/blob/main/samples/sample.ipynb

## 参考資料

位置情報はいろいろ難しいところがあり、簡単と思っていても処理を失敗することもあるので、ドメイン知識が重要です。私も時折やらかすため、いろいろ勉強中です。

最後に今回の記事を読むうえで役立つと思う資料をピックアップしておきます。

座標参照系に関して知りたい方はQGISのドキュメントをご参照ください

https://docs.qgis.org/3.28/ja/docs/gentle_gis_introduction/coordinate_reference_systems.html

平面直角座標系に関しては国土地理院の解説ページをご参照ください

https://www.gsi.go.jp/sokuchikijun/jpc.html