ユーザーがかなり限定されますが...
※ テストが十分でないのでエラー等指摘いただけると助かります。
lemファイルとは
- lemファイルとは、航空LiDARの数値標高データの形式の一つです。
- lemファイルとその属性情報を示すcsvファイルがセット。
- フォーマットは次を参照
https://gigaplus.makeshop.jp/jmcnet/download/Form2m.txt
pythonスクリプト
- lem2geotiff.pyの関数translateが本体です。
処理の流れは以下のとおりです。
- csvファイルより必要情報を読み取る。
- lemファイルを読み、ndarrayに変換する。
- ndarrayからxarrayを作成する。
- 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ファイルを使うことが多い。
- pythonの場合
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
参考サイト
raster - Python equivalent of gdalbuildvrt - Geographic Information Systems Stack Exchange
