実行環境
Windows10
Python 3.9.5
geopandas 0.10.2
pandas 1.3.1
shapely 1.8.0
はじめに
水理学のいろはも知らなければ iRIC も未経験なのですが、ちょっと機会があって iRIC格子CSVファイルを GIS データに変換するということをしました。
今回はそのとき python で書いたコードを載せたいと思います。
iRIC は、CGNS ファイルや VTK ファイルなどの形式で格子をエクスポートすることができます。
ファイルフォーマットについてはこちらに載っています。
これらのファイルを扱うライブラリで GIS データに変換できないだろうかと少し調べましたが、ちょっと分かりませんでした。
ググっても iRIC の格子を一般的な GIS データに直接変換する方法が見つかりませんでした。
きっとスマートなやり方があるんだろうとは思いますが、今回は愚直に変換していきます。
格子について
個人的にはメッシュという呼び方の方がなじみがあるのですが、iRIC のマニュアルに準じてここでも格子と表記します。
iRIC格子CSVファイルでは BFC 格子を扱います。
流れの方向に沿って並んだ4角形の格子になります。
仕様はこちらをご確認ください。
I が上流から下流に向かう方向、J が右岸から左岸へ向かう方向のようです。
K という次元も扱えるようですが、IJ の2次元で考えます。
また、Nays2DH というソルバーの格子データを想定しています。
他のソルバーでも使えるとは思いますが、テストはしてないです。
サンプルファイルについて
こちらのチュートリアルに従って私が作成した格子を使います。
サンプルファイル（1901行のうち先頭10行）
IMAX,JMAX,KMAX
73,26,1
I,J,K,X,Y,Z,N_Elevation,N_Elevation_zb,C_Obstacle,C_Fix_movable,C_vege_density,C_vege_height,C_roughness_cell,C_mix_cell
0,0,0,-5.18970336e+04,-9.41273264e+04,0,1.01386906e+01,0.00000000e+00,0,1,0.00000000e+00,0.00000000e+00,4.00000000e-02,0
1,0,0,-5.18979442e+04,-9.42171220e+04,0,9.61127733e+00,0.00000000e+00,0,1,0.00000000e+00,0.00000000e+00,4.00000000e-02,0
2,0,0,-5.19007733e+04,-9.43088529e+04,0,9.20632464e+00,0.00000000e+00,0,1,0.00000000e+00,0.00000000e+00,4.00000000e-02,0
3,0,0,-5.19056737e+04,-9.44015579e+04,0,8.67082308e+00,0.00000000e+00,0,1,0.00000000e+00,0.00000000e+00,4.00000000e-02,0
4,0,0,-5.19131691e+04,-9.44945179e+04,0,7.97978274e+00,0.00000000e+00,0,1,0.00000000e+00,0.00000000e+00,4.00000000e-02,0
5,0,0,-5.19242780e+04,-9.45871201e+04,0,7.22250258e+00,0.00000000e+00,0,1,0.00000000e+00,0.00000000e+00,4.00000000e-02,0
6,0,0,-5.19404711e+04,-9.46787950e+04,0,6.47999352e+00,0.00000000e+00,0,1,0.00000000e+00,0.00000000e+00,4.00000000e-02,0
...
ソルバーマニュアル iRIC v3対応
北海道の石狩川です。
iRIC で表示しています。（背景は地理院地図の淡色地図）
上図では右上が上流になります。
格子の形状は水理計算の精度や安定性に関わる重要な要素ではあると思いますが、今回は GIS データへの変換が目的なので、多少いびつなところがあってもお許しください。
出力について
格子点をポイントデータ、セルをポリゴンデータで出力します。
また、格子点の属性値をポイントデータに、セルの属性値をポリゴンデータに付与します。
仕様にあるとおり、属性名の先頭が "N_" なら格子点、 "C_" ならセルの属性です。
セルのデータ数は、格子点のデータ数に対して IJ 各列の終端が削られる分少なくなります。
ファイル形式はシェープファイルを標準としていますが、geopandas の GeoDataFrame.to_file() を使っているので、GeoJSON などほかの GIS データファイル形式を指定することもできます。
コード
import os.path
import geopandas as gpd
import pandas as pd
from shapely.geometry import Polygon
class IricMesh:
def __init__(self, fp: str, epsg=None, colz='N_Elevation') -> None:
""" iRIC格子CSVファイルのデータ
Parameters
--------
fp : str
iRIC格子CSVファイルのパス
epsg : int
EPSGコード
デフォルトは None で座標系設定なし
colz : str
格子点の高さを格納した列の名前
デフォルトは 'N_Elevation'
Attributes
--------
imax : int
I方向（縦断方向）格子点数
jmax : int
J方向（横断方向）格子点数
kmax : int
K方向格子点数（使用しない）
df : pd.DataFrame
CSVファイルの情報とジオメトリを格納
epsg : int
EPSGコード
colz : str
格子点の高さを格納した列の名前
"""
with open(fp, 'r') as f:
f.readline()
ijk = [int(n) for n in f.readline().strip().split(',')]
self.imax, self.jmax, self.kmax = ijk
self.df = pd.read_csv(fp, skiprows=2)
self.df = self.df.sort_values(by=['J', 'I']) # 念のためソート
self.df[['point', 'polygon']] = None # 格子点とセルのジオメトリ列
self.epsg = epsg
self.colz = colz
def point(self) -> None:
""" 格子点のジオメトリを生成
"""
self.df['point'] = gpd.points_from_xy(
self.df['X'], self.df['Y'], self.df[self.colz])
def polygon(self) -> None:
""" セルのジオメトリを生成
"""
df = self.df
# 格子点以外の残りの３つの頂点を列として追加する
df[['x2', 'y2', 'z2']] \
= df.loc[:, ['X', 'Y', self.colz]].shift(-1)
df[['x3', 'y3', 'z3']] \
= df.loc[:, ['X', 'Y', self.colz]].shift(-self.imax-1)
df[['x4', 'y4', 'z4']] \
= df.loc[:, ['X', 'Y', self.colz]].shift(-self.imax)
def func_polygon(row):
""" ４つの頂点から生成されるポリゴンのジオメトリを返す関数
"""
x1, y1, z1 = row['X'], row['Y'], row[self.colz]
x2, y2, z2 = row['x2'], row['y2'], row['z2']
x3, y3, z3 = row['x3'], row['y3'], row['z3']
x4, y4, z4 = row['x4'], row['y4'], row['z4']
polygon = Polygon(
[(x1, y1, z1), (x2, y2, z2), (x3, y3, z3), (x4, y4, z4)]
)
return polygon
# IJ各列の終点にはセルがないので除外して、ジオメトリを生成
df.loc[
(df['I'] != self.imax - 1) & (df['J'] != self.jmax - 1), 'polygon'
] = df.apply(func_polygon, axis=1)
# 不要な列の削除
df = df.drop(
['x2', 'y2', 'z2', 'x3', 'y3', 'z3', 'x4', 'y4', 'z4'],
axis=1
)
self.df = df
def output(self, fp_point: str, fp_polygon: str, driver=None) -> None:
""" GISファイルへの出力
Parameters
--------
fp_point : str
格子点の出力ファイルパス
fp_polygon : str
セルの出力ファイルパス
driver : str
出力フォーマット
デフォルトの None だとシェープファイル形式
GeoJSON 形式なら 'GeoJSON' を指定
"""
def geom(prefix: str, colgeom: str) -> gpd.GeoDataFrame:
""" ファイル出力用のジオデータフレームを返す関数
Parameters
--------
prefix : str
属性名の接頭辞
格子点の属性なら"N_"、セルの属性なら"C_"
colgeom : str
ジオメトリを格納した列の名前
"""
colall = self.df.columns.to_list()
# 必要な属性の名前のリスト
cols = ['I', 'J', 'K', 'X', 'Y', 'Z'] \
+ [col for col in colall if col.startswith(prefix)] \
+ [colgeom]
gdf = gpd.GeoDataFrame(
self.df[cols], geometry=colgeom, crs=self.epsg)
return gdf
geom(prefix='N_', colgeom='point').to_file(fp_point, driver=driver)
geom(prefix='C_', colgeom='polygon').to_file(fp_polygon, driver=driver)
if __name__ == '__main__':
fold = os.path.dirname(__file__)
fp_csv = os.path.join(fold, 'mesh.csv')
fp_pt = os.path.join(fold, 'shp', 'mesh_point.shp')
fp_poly = os.path.join(fold, 'shp', 'mesh_polygon.shp')
code = 2454 # JGD2000 / Japan Plane Rectangular CS XII
iricmesh = IricMesh(fp=fp_csv, epsg=code)
iricmesh.point()
iricmesh.polygon()
iricmesh.output(fp_point=fp_pt, fp_polygon=fp_poly)
フィールド名が10文字を超えるので、実行時 UserWarning が出ます。
GISソフトでチェック
ArcGIS Pro で出力ファイルを表示してみます。
平面表示の背景は地理院地図の淡色地図です。
立体表示では高さ方向を5倍に強調し、格子点を高さで色付けしています。
平面表示
立体表示（下流から上流を見る）
ポイントデータ属性表示
ポリゴンデータ属性表示
シェープファイルの制限でフィールド名が10文字になります。
また、出力時に型指定をしていないのですべてのフィールドが Double 型です。
形状や座標については問題なさそうです。
さいごに
iRIC格子CSVファイルを GIS データに変換するコードをご紹介しました。
しかしながら、やはりほかにもっと一般的な方法があると思われます。