趣味で計算流砂水理

趣味で計算流砂水理 Computational Sediment Hydraulics for Fun Learning

数値計算とか河川工学とかプログラミングのことを書いています

MENU

オープンデータのALB測量よりjupyter上で陰影図を作成する

スポンサーリンク

こんなのを作りました。

f:id:SedimentHydraulics:20210615104742g:plain

モチベーション

グリーンレーザありの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

陰影図を作成

陰影図は以下のサイトを参照して作成しました。

desktop.arcgis.com

変換式は次のとおりです。

\begin{align}
slope &= \arctan{ \left( \sqrt{ \left( \frac{\partial z}{\partial x} \right)^2 + \left( \frac{\partial z}{\partial y} \right)^2}\right) } \\
aspect &= \arctan2{  \left( \frac{\partial z}{\partial y} , \frac{\partial z}{\partial x} \right) } \\
zenith &= \pi/2 -  altitude \\
azimuth &= 2\pi -  azimuth + \pi/2 \\
hillshade &= 255  \left\{ \cos( zenith) \cos(slope) +  \sin(zenith) \sin( slope) \cos(azimuth - aspect) \right\} 
\end{align}

ここに、slope:最急勾配、aspect:最急勾配の方位角、altitude:光源の高度、azimuth:光源の方位角

高速化の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

github.com


gif作成

computational-sediment-hyd.hatenablog.jp