### 平面直角座標系を求めるツールって意外となくない？
---
平面直角座標系とは
---
---
日本を19のゾーンに分割し、各ゾーンに座標原点を設けて平面として捉えた座標系
---
系番号と座標を指定して緯度経度に変換するツール
https://www.n-survey.com/online/gmap2.htm
---
---
---
そもそも指定した緯度経度がどこの座標系なのか理解していないといけない。
---
- 9系/360613.58925/1400516.27815
---
---
- 16系/360613.58925/1400516.27815
---
---
座標系を判別するプログラムはなさそう？
---
ないなら作ればいいじゃない
---
https://github.com/nokonoko1203/search_zone_number
---
市区町村名・緯度経度から平面直角座標系を取得できます。
---
- 緯度経度で検索
```bash=
% pipenv run python -m search_zone_number lnglat -n 134.814283 -a 35.543052
GST_CSS_NAME system_number
1676 豊岡市 5
```
---
- 市区町村名で検索
```bash=
% pipenv run python -m search_zone_number name -n 豊岡市
GST_CSS_NAME system_number
1676 豊岡市 5
```
---
実装
---
- まずはデータを作成
- pref_codeに都道府県コードを入力すると境界データが取得できる
```
https://www.e-stat.go.jp/gis/statmap-search/data?dlserveyId=A002005212015&code={pref_code}&coordSys=1&format=shape&downloadType=5
```
---
- shp -> gdf -> gpkg
- gpd.read_file()はzipファイルを直で読める
- 地物の融合も簡単
```python=
zip_path_list = [Path(l) for l in glob.glob("./*.zip")]
town_gdf = pd.concat([
gpd.read_file("zip://" + str(zipfile))
for zipfile in zip_path_list
]).pipe(gpd.GeoDataFrame)
town_gdf["AREA_CODE"] = town_gdf['PREF'].str.cat(town_gdf['CITY'])
city_gdf = town_gdf.dissolve(by="AREA_CODE", as_index=False)
pref_gdf = city_gdf.dissolve(by="PREF", as_index=False)
pref_gdf.to_file("boundary.gpkg", layer='pref', driver="GPKG")
city_gdf.to_file("boundary.gpkg", layer='city', driver="GPKG")
town_gdf.to_file("boundary.gpkg", layer='town', driver="GPKG")
```
---
- 全て手作業で系番号を設定…
```python=
# 1系
town_gdf['system_number'] = None
one_city_list = [
"奄美市",
"名瀬市",
"十島村",
...
town_gdf["system_number"].mask((town_gdf["CITY_NAME"].isin(one_city_list)) & (
town_gdf["PREF_NAME"] == "鹿児島県"), 1, inplace=True)
town_gdf["system_number"].mask(town_gdf["PREF_NAME"] == "長崎県", 1, inplace=True)
```
---
- ファイルに書き出し
```python=
import geopandas as gpd
gdf = gpd.read_file("merge_boundary.gpkg", layer="add_system_number")
gdf["GST_CSS_NAME"] = gdf['GST_NAME'].str.cat(gdf['CSS_NAME'], na_rep="")
dissolve_key = ["GST_CSS_NAME", "system_number"]
city_boundary = gdf.dissolve(by=dissolve_key, as_index=False)
city_boundary.to_feather("merge_city_boundary.feather")
city_boundary.to_file("merge_city_boundary.gpkg", layer="merge_city_boundary", driver="GPKG")
```
---
- featherファイル？
- shp/csvなどより高速に読み書きできる形式
- 最近GeoPandasに取り込まれた
- gpd.readfile()で読み込める
- [https://geopandas.readthedocs.io/en/latest/docs/user_guide/io.html](https://geopandas.readthedocs.io/en/latest/docs/user_guide/io.html)
---
- 町丁目レベル
---
- わかりやすいように１９座標系にマージ
---
- なんやかんやで市区町村レベルに結合（これを利用！）
---
- こんな感じのデータを保持
---
データがあればちょちょいのちょい
---
- 市区町村名で検索
```python=
gdf[gdf["GST_CSS_NAME"].str.contains(city_name)]
```
---
- 緯度経度で検索
```python=
point = gpd.GeoSeries([Point([(lng, lat)])])
point.set_crs(epsg=4612, inplace=True)
gdf[gdf["geometry"].apply(lambda x: point.within(x).any()) == True]
```
---
- 今後
- WebAPIにしたい
- PyPIに登録しときたい
---
完
{"metaMigratedAt":"2023-06-15T18:32:22.628Z","metaMigratedFrom":"Content","title":"Untitled","breaks":true,"contributors":"[{\"id\":\"501938f7-150f-4a68-b27a-4907062f3c90\",\"add\":3670,\"del\":129}]"}