実行環境
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行）

mesh.csv

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 データに変換するコードをご紹介しました。

しかしながら、やはりほかにもっと一般的な方法があると思われます。