趣味で計算流砂水理

趣味で計算流砂水理

Computational Sediment Hydraulics for Fun

MENU

地理情報付きラスタ画像から任意座標のデータを抽出するpythonスクリプト

スポンサーリンク

最近よく使うスクリプトなので備忘録的にまとめました。


この記事のポイント
  • ラスタ画像から任意座標のデータの抽出はrioxarrayを使うと超簡単
  • この手の方法は航空機、ドローン等の高密度測量データから計算データセットを作成する際には必須

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:])
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を使用する。

データの取得は以下のライブラリを使用する。

libgsidem2el

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

GitHub - computational-sediment-hyd/ExtractingPointDataFromRasterImage: Extracting A Point Data From A Raster Image


  • グラフの作成はHoloviewsを使用

computational-sediment-hyd.hatenablog.jp

  • はてブロへのhtmlグラフの埋め込みは以下を参照

computational-sediment-hyd.hatenablog.jp