最近よく使うスクリプトなので備忘録的にまとめました。
この記事のポイント
- ラスタ画像から任意座標のデータの抽出はrioxarrayを使うと超簡単
- この手の方法は航空機、ドローン等の高密度測量データから計算データセットを作成する際には必須
- environment
- サンプルデータ
- テキストデータからgeotiffを作成
- 任意座標のデータをgeotiffから抽出する
- 応用編：左右岸杭座標を入力してgeotiffから河川横断データを作成する
- おまけ：国土地理院の標高データと比較
- Github
environment
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
サンプルデータ
那賀川(code:8808060001)のALBデータより 04GG542_0.5gを使用する。
テキストデータからgeotiffを作成
過去記事(欠損したグリッド標高テキストデータからgeotiffを作成するpythonスクリプト - 趣味で計算流砂水理)と同様。
データサイズを小さくするため、単位系をmからcmに変換し、geotiffをint32型とする変更を加えた。
作成したファイルは「04gg542_0.5g.tif」
こちらからダウンロード可能。
任意座標のデータをgeotiffから抽出する
import numpy as np import xarray as xr
- 任意の点の平面直角座標xy
pl = [98891.63181796, 104393.37756069]
- ラスタの情報を読み込む
# read raster ds = xr.open_rasterio('04gg542_0.5g.tif') # get raster properties dy, dx = ds.res xOrigin =ds.x.values[0] - 0.5*dx yOrigin =ds.y.values[0] + 0.5*dy ynum = ds.shape[1] xnum = ds.shape[2] z = ds.values[0] nepsg = int(ds.attrs['crs'][-4:])
- 任意の点のgeotiff上のindexを取得 : 過去記事(欠損したグリッド標高テキストデータからgeotiffを作成するpythonスクリプト - 趣味で計算流砂水理)のラスタ作成の逆の処理
xind = np.floor( (pl[0] - xOrigin)/dx ).astype(int) yind = np.floor( (yOrigin - pl[1])/dy ).astype(int) print(z[yind, xind]/100)
21.7
- 図で確認:web表示の関係でラスタの画素が粗くなっているがローカルでは0.5m解像度になる。
作成したコードはこちら
応用編：左右岸杭座標を入力してgeotiffから河川横断データを作成する
- 左右岸の杭座標（平面直角座標系）
# 左岸杭座標 pl = [98891.63181796, 104393.37756069] # 右岸杭座標 pr = [ 99438.53777557, 104214.89229045]
- 測線上にdelta間隔の点群pointsを発生させる。
- 杭の外側（堤内地側）にLex延伸する。
pl = np.array(pl) pr = np.array(pr) # 測線の長さ L = np.linalg.norm(pr - pl) # 単位ベクトルe e = (pr - pl)/ L # 杭の外側（堤内地側）にLex延伸する。 Lex = 50 # 測線上にdelta間隔の点群を発生させる。 delta = 0.5 dl = np.arange(-Lex, L+Lex, delta) points = np.array( [pl+e*dlp for dlp in dl] )
- points のgeotiff上のindexを取得
indx = np.floor( (points[:,0] - xOrigin)/dx ).astype(int) indy = np.floor( (yOrigin - points[:,1])/dy ).astype(int) Z = np.array([z[iy, ix]/100 for ix,iy in zip(indx, indy)])
- 図化してチェック
import holoviews as hv hv.extension('bokeh')
g = hv.Curve((dl,Z))
hv.save(g,'ALBcrosssect.html')
g
- 平面図で確認:web表示の関係でラスタの画素が粗くなっているがローカルでは0.5m解像度になる。
作成したコードはこちら
おまけ：国土地理院の標高データと比較
国土地理院の標高データのうちは5m間隔のLPデータから作成されているDEM5Aを使用する。
データの取得は以下のライブラリを使用する。
import libgsidem2el as gsi import geopandas as gpd from shapely.geometry import Point
- 測線上にdelta間隔の点群を発生させる。
delta = float(5) dl2 = np.arange(-Lex, L+Lex, delta) points2 = np.array( [pl+e*dlp for dlp in dl2] )
- 国土地理院のデータは緯度経度で検索する必要があるため、座標変換を行なう。
- geopandasを使う方法を示すが他にもいろいろな方法がある。
gdf = gpd.GeoDataFrame(geometry=[Point(p[0],p[1]) for p in points2], crs={'init': 'epsg:'+str(nepsg)}) # geopandas : change coordinate gdf = gdf.to_crs(epsg=6668)
- 上記のライブラリを使って標高を取得
dem5 = gsi.libgsidem2el('DEM5A') Z2 = [] for go in gdf.geometry.values: zt = dem5.getEL(go.x, go.y, zoom = 15) if(zt == 'e') or (zt == 'outside') : zt = float(-9999) Z2.append(zt) Z2 = np.array(Z2, dtype=float)
- 図化してチェック
ALB:0.5mメッシュとLP:5mメッシュで大局的にはそんなに変わらない気がします。
河岸とか堤防とかをみる場合は細かい方がいいかもという感じです。
g = hv.Curve((dl,Z), label='ALB') * hv.Curve((dl2,Z2), label='LP') hv.save(g,'ALBandLP.html') g
Github
- グラフの作成はHoloviewsを使用
computational-sediment-hyd.hatenablog.jp
- はてブロへのhtmlグラフの埋め込みは以下を参照
