とてつもなくユーザーが少ない話ですが。。。
目的
航空測量成果に納められているグリッド標高テキストデータは、 データ容量節約のため、任意の図郭内で測量対象範囲外を除いたデータのみが記載されていることが多いです。 この場合、データをリシェープするだけではラスタ形式（二次元配列）を作ることが出来ないので工夫が必要です。
なお、lemファイルがある場合は次の記事を参照下さい。
computational-sediment-hyd.hatenablog.jp
ファイル名から図郭北西座標を取得
参考サイト :国土基本図図郭とは｜図郭コードの読み方 - 空間情報クラブ｜インフォマティクス運営のメディアサイトをみると、 ファイル名から図郭を取得可能なんですが、測量精度等により結構個性が出たりします。
例えば、今回サンプルとして扱ったデータの図郭名はこんな感じです。
そのため、汎用化は難しいですが、今回の例に対応した図郭北西座標を求めるスクリプトは次のようになります。
特に複雑な処理は無いですがord関数で文字列を数値化している点が特徴でしょうか。
ファイル名が小文字の場合は、ラムダ式内をord('a')に変更します。
ファイル名が小文字の場合は、isCapsオプションをFalseに変更します。
import re def getcoords(word, xorg, yorg, dx, dy, isCaps=True): # anum = lambda x : ord(x) - ord('A') if isCaps: anum = lambda x : ord(x) - ord('A') else: anum = lambda x : ord(x) - ord('a') v = str(word) if re.match('\d+', v): iy, ix = int(v[0]), int(v[1]) else: iy, ix = anum(v[0]), anum(v[1]) yc = yorg - iy*dy xc = xorg + ix*dx return xc, yc # 4分割 def getcoords2(word, xorg, yorg, dx, dy): v = int(word) if v == 1: iy, ix = 0, 0 elif v == 2: iy, ix = 0, 1 elif v == 3: iy, ix = 1, 0 elif v == 4: iy, ix = 1, 1 else: print('error') yc = yorg - iy*dy xc = xorg + ix*dx return xc, yc
w = '09LD3312' nepsg = int(w[:2]) # 図郭の北西の座標を求める # 地図情報レベル 50000 xc, yc = getcoords(w[2:4], xorg=-160000, yorg=300000, dx=40000, dy=30000) # 地図情報レベル 5000 xc, yc = getcoords(w[4:6], xorg=xc, yorg=yc, dx=4000, dy=3000) # 地図情報レベル 5000を4分割 xc, yc = getcoords2(w[6], xorg=xc, yorg=yc, dx=2000, dy=1500) # さらに4分割 xc, yc = getcoords2(w[7], xorg=xc, yorg=yc, dx=1000, dy=750)
元データをラスタ用の二次元配列に変更
上記の図郭北西座標をメッシュサイズよりXY座標のインデックスを取得し、標高Zの二次元配列に与える。 今回の例ではメッシュサイズは1m、欠測値は-9999としています。
import numpy as np import pandas as pd df = pd.read_csv('09LD3312_1g.txt', header=None) X = df[1].values Y = df[2].values Z = df[3].values delta = float(1.0) dx=1000 dy=750 ix = ((X - xc - 0.5*delta)/delta).astype(int) iy = ((- Y + yc - 0.5*delta)/delta).astype(int) zout= np.full((int(dx/delta),int(dy/delta)), float(-9999) ) for ixp, iyp, zp in zip(ix, iy, Z) : zout[ixp,iyp] = zp
データをgeotiff化
ここからのこのサイトでは定番の流れでxarray → rioxarrayとしてgeotiffでエクスポートします。 ポイントは出力する二次元配列を（y,x）の順に定義することです。 複数のバンドを持つ場合は、coordsにbandを加えて（band,y,x）の三次元配列にすること。
import xarray as xr import rioxarray xcoord = np.arange(xc, xc+dx, delta) + 0.5*delta ycoord = np.arange(yc, yc-dy, -delta) - 0.5*delta epsg = str(6668 + nepsg) ds = xr.Dataset({'z': (['y','x'], zout.T) }, coords={'x': xcoord, 'y': ycoord}) ds = ds.rio.write_crs('EPSG:' + epsg, inplace = True) # export geotiff file out = ds['z'].rio.to_raster('09LD3312_1g.tif') del out
以上で完成です。
おまけ：VRTファイルの作成
geotiffは膨大なファイル数になることが多いので、ソフトウェアで取り扱う場合はVRTファイルを使うことが望ましい。
from osgeo import gdal # gdal 2.1以降 my_vrt = gdal.BuildVRT('output.vrt', glob.glob( '*.tif'), VRTNodata=-9999, srcNodata=-9999) my_vrt = None
参照：Pythonから使いやすくなったGDAL 2.1がステキだ | giscience.jp
Environment
rasterとxarrayの変換にはrioxarrayモジュールを使用している。
詳細を省略するが、gdalを使用するよりもコード量が大幅に少なくなるため、非常に便利である。
そのための仮想環境tmpの構築はanacondaを用いると次のようになる。
conda create -y -n tmp python=3.7 conda activate tmp conda install -y -c conda-forge rioxarray conda install -y -c conda-forge gdal conda install -y -c conda-forge jupyter
