こんなのを作りました。
モチベーション
グリーンレーザありのALB(Airborne Laser Bathymetry)に河川測量成果について、 九頭竜川(code:8606070001) 、那賀川(code:8808060001) でオープンデータ化されて、昨年(一昨年?)に少しだけ話題になりました。
グリーンレーザの測量精度については思うところはありますが本題とはずれるのでまたの機会に。
この手のデータを可視化する場合、陰影図をよく使いますが、GISソフトを開くのが面倒なのでjupyter上で書いてみようかなと。
使用データ
上記の那賀川のデータより、50cmメッシュデータの「04gg542_0.5g」図郭のデータを使いました。
グリットテキストデータよりgeotiffを作成
この方法は、以下の記事に詳述しているため、結果のコードのみ記載します。 04gg542_0.5g.txtから04gg542_0.5g.tifを作りました。
computational-sediment-hyd.hatenablog.jp
import re import glob import pandas as pd import numpy as np import xarray as xr import rioxarray import os
def getcoords(word, xorg, yorg, dx, dy): 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
fs =glob.glob('../data/*.txt') for fp in fs: f = os.path.basename(fp) df = pd.read_csv(fp, header=None) X = df[1].values Y = df[2].values Z = df[3].values w = f[:8] 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) delta = float(0.5) dx=2000 dy=1500 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 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) # ds.rio.reproject("epsg:****") # export geotiff file out = ds['z'].rio.to_raster( '../' + f[:-4] + '.tif') del out
陰影図を作成
陰影図は以下のサイトを参照して作成しました。
変換式は次のとおりです。
ここに、:最急勾配、:最急勾配の方位角、:光源の高度、:光源の方位角
高速化のnumbaを使っています。
computational-sediment-hyd.hatenablog.jp
import numpy as np import xarray as xr import hvplot.xarray from numba import jit
@jit(nopython=True) def Hillshade(z, dy, dx, Altitude = 45.0 * np.pi/180, Azimuth = 315.0 * np.pi/180): # Altitude : 光源高度 # Azimuth : 光源の方位角 Zenith = np.pi/2 - Altitude #天頂角 Azimuth = 2*np.pi - Azimuth + np.pi/2 if Azimuth >= 2*np.pi : Azimuth - 2*np.pi z_factor = 1.0 # zの強調 ny, nx = z.shape Hillshade = np.full_like(z, 0, dtype=np.uint8) for j in range(1, ny-1): for i in range(1, nx-1): if np.isnan(z[j,i]) : tmp = int(0) else: dzdx = 0.25*( 0.5*(z[j+1,i+1] - z[j+1,i-1])/dx + 1.0*(z[j ,i+1] - z[j ,i-1])/dx + 0.5*(z[j-1,i+1] - z[j-1,i-1])/dx) dzdy = 0.25*( 0.5*(z[j+1,i+1] - z[j-1,i+1])/dy + 1.0*(z[j+1,i ] - z[j-1,i ])/dy + 0.5*(z[j+1,i-1] - z[j-1,i-1])/dy) Slope = np.arctan( z_factor * np.sqrt(dzdx**2 + dzdy**2) ) if dzdx == 0: if dzdy > 0 : Aspect = np.pi / 2.0 elif dzdy < 0 : Aspect = 2.0 * np.pi - np.pi / 2.0 else : Aspect = 0.0 else : Aspect = np.arctan2(dzdy, -dzdx) if Aspect < 0 : Aspect = 2.0 * np.pi + Aspect # Hillshade tmp = 255.0 * ((np.cos(Zenith) * np.cos(Slope)) + (np.sin(Zenith) * np.sin(Slope) * np.cos(Azimuth - Aspect))) tmp = int(tmp) Hillshade[j,i] = tmp return Hillshade
dstif = xr.open_rasterio('../04gg542_0.5g.tif') z = dstif.values z = np.where(z==-9999, np.nan, z) z = z[0] dy, dx = dstif.attrs['res'] a = Hillshade(z, dy, dx)
図化
定番のpyvizを使った方法です。
computational-sediment-hyd.hatenablog.jp
import numpy as np import xarray as xr import cartopy.crs as crs import hvplot.xarray import geoviews as gv from numba import jit gv.extension('bokeh', logo=False)
dd = xr.DataArray(a , coords={'x':dstif['x'], 'y':dstif['y']} , dims=['y', 'x'] , attrs= dstif.attrs ) nepsg = int(dstif.attrs['crs'][-4:]) # html export # g = dd.hvplot.image(x='x',y='y', crs=crs.epsg(nepsg), cmap='gray',geo=True # , dynamic=False, rasterize=False,clim=(1,255), colorbar=False, alpha=0.7, frame_width=600).options(clipping_colors={'min':'transparent'}) # on jupyter g = dd.hvplot.image(x='x',y='y', crs=crs.epsg(nepsg), cmap='gray',geo=True , dynamic=True, rasterize=True, clim=(1,255), colorbar=False, alpha=0.7, frame_width=600).options(clipping_colors={'min':'transparent'}) back = gv.WMTS('https://mt1.google.com/vt/lyrs=s&x={X}&y={Y}&z={Z}', name="GoogleMapsImagery", crs=crs.epsg(6672)) gall = back * g
jupyter版とhtml出力版がありますが、後者ではブログ表示用にhtmlにデータを埋め込むため、図化方法を静的にしています。 jupyterで使う場合は動的のほうが軽量のため扱いやすいです。
d = hvplot.save(gall,'hillshade.html') del d
github
gif作成