趣味で計算流砂水理

Computational Sediment Hydraulics for Fun

lemファイルをgeotiffファイルに変換するpythonスクリプト

ユーザーがかなり限定されますが...

※ テストが十分でないのでエラー等指摘いただけると助かります。


lemファイルとは

  • lemファイルとは、航空LiDARの数値標高データの形式の一つです。

参照:https://www.harrisgeospatial.co.jp/Portals/74/VIS_Japan/documents/3_ENVI2017_CaseStudy_Toda.pdf?ver=2017-11-09-141636-000#page=24

  • lemファイルとその属性情報を示すcsvファイルがセット。
  • フォーマットは次を参照

https://gigaplus.makeshop.jp/jmcnet/download/Form2m.txt

pythonスクリプト

  • lem2geotiff.pyの関数translateが本体です。
  • 処理の流れは以下のとおりです。

    1. csvファイルより必要情報を読み取る。
    2. lemファイルを読み、ndarrayに変換する。
    3. ndarrayからxarrayを作成する。
    4. xarrayをgeotiffとして出力
  • 引数は次の4つのみ。filepath以外はオプショナルです。

    • filepath : str - lemファイルのフルパス
    • valm1111 : float, default 60000 - データ作成範囲外の場合のフラグ-1111の値
    • valm9999 : float, default 50000 - 海部及び陸水部の場合のフラグ-9999の値
    • dtype : type, default np.uint16 - geotiffファイルのデータ型
def translate(filepath, valm1111=60000, valm9999=50000, dtype=np.uint16 ):
    fpath = Path(filepath)
    
    # reading csv file
    df = pd.read_csv(fpath.with_suffix('.csv'), encoding='ShiftJIS',header=None, index_col=0)
    
    nx = int(df.loc['東西方向の点数'].values[0])
    ny = int(df.loc['南北方向の点数'].values[0])
    y1 = float(df.loc['区画左下X座標'].values[0])
    x0 = float(df.loc['区画左下Y座標'].values[0])
    y0 = float(df.loc['区画右上X座標'].values[0])
    x1 = float(df.loc['区画右上Y座標'].values[0])
    dx = float(df.loc['東西方向のデータ間隔'].values[0])
    dy = float(df.loc['南北方向のデータ間隔'].values[0])
    nepsg = int(df.loc['平面直角座標系番号'].values[0])
    
    # reading lem file
    with open(filepath, "r") as f:
        lines = f.readlines()
    
    # Important: The dimensions of the array must be in the order from y to x.
    arr = np.empty((ny, nx))
    for i, l in enumerate(lines) :
        ll = l.replace('\n','')
        s = ll[10:]
        v = [s[i: i+5] for i in range(0, len(s), 5)]
        arr[i,:] = np.array(v)
    
    # make xarray
    val = np.array(arr, np.float)
    
    # Changing the flag -9999 for the sea and land and water sections and the flag -1111 outside the data creation area into any given value.
    # dafault : -9999 to 50000, -1111 to 60000
    val = np.where(val == -1111, valm1111, val)
    val = np.where(val == -9999, valm9999, val)
    
    # changing data type(default uint16)
    val = val.astype(dtype)
    
    dx *= 100
    dy *= 100
    xarr = np.arange(x0,x1,dx)
    xarr = xarr + dx/2
    yarr = np.arange(y0,y1,-dy)
    yarr = yarr - dy/2
    
    epsg = str(6668 + nepsg)
    
    ds = xr.Dataset({'z': (['y','x'], val) }, coords={'x': xarr/100, 'y': yarr/100}) 
    ds = ds.rio.write_crs('EPSG:' + epsg, inplace = True)
    # ds.rio.reproject("epsg:****")
    
    # export geotiff file
    ds['z'].rio.to_raster( fpath.stem + '.tif')

Usage

Exmaple 1

単一ファイルの場合

import lem2geotiff as l2g

l2g.translate('lem/****.lem')

Exmaple 2

複数ファイルの一括処理

import lem2geotiff as l2g
import glob

for ff in glob.glob( 'lem/*.lem'):
    l2g.translate(ff)
おまけ

lemファイルは膨大なファイル数になることが多いので、ソフトウェアで取り扱う場合はVRTファイルを使うことが多い。

from osgeo import gdal 

my_vrt = gdal.BuildVRT('output.vrt', glob.glob( '*.tif'))
my_vrt = None
  • gdalの場合
gdalbuildvrt out.vrt *.tif

Environment

rasterとxarrayの変換にはrioxarrayモジュールを使用している。

詳細を省略するが、gdalを使用するよりもコード量が大幅に少なくなるため、非常に便利である。

Anaconda

仮想環境lemを作成する。

conda create -y -n lem python=3.7
conda activate lem 
conda install -y -c conda-forge rioxarray
conda install -y -c conda-forge gdal
conda install -y -c conda-forge jupyter

Google Colaboratory

  • rioxarrayをインストールする。
  • ドライブをマウントして作業ディレクトリに移動する。
!pip install rioxarray
from google.colab import drive
drive.mount('/content/drive')
%cd "/content/drive/My Drive/*"

github

github.com

参考サイト

raster - Python equivalent of gdalbuildvrt - Geographic Information Systems Stack Exchange

[Python]可読性を上げるための、docstringの書き方を学ぶ(NumPyスタイル) - Qiita

docstringのstyle3種の例 - Qiita

numpydoc docstring guide — numpydoc v1.2.dev0 Manual