ミャンマー某河川の河口部：ETOPO1 + geoviews - 趣味で計算流砂水理 Computational Sediment Hydraulics for Fun Learning

またまた地図ネタです．

先日ネタにでましたミャンマーのある河川の河口部について，ETOPO1とgeoviewsで絵を書いてみました．

大きいものはこちら

スケールがわかりにくいですが，河口から300km程度は，10m程度の浅い水深の範囲が続いています．
河口から50kmくらいがぐちゃぐちゃな感じです．この辺りが活発に河床変動が生じているような気がします．
スケールが大きくて面白いです．
landsatのデータもかなり面白いので，後日絵を書いてみます．

ETOPO1 + geoviewsの作り方

驚くほど簡単でした．手順をまとめておきます．

モジュールのロード

以下のとおり，geoviewsとxarrayを読みます．

xarrayは多次元データの解析用のモジュールで非常に高速です．後日記事を書きます．

import xarray as xr
import geoviews as gv
import geoviews.tile_sources as gts

gv.extension('bokeh')

ETOPO1をロード

以下から，ETOPO1のnetCDFデータをダウンロードします．これが一番時間がかかります．

https://www.ngdc.noaa.gov/mgg/global/relief/ETOPO1/data/ice_surface/grid_registered/netcdf/

リードは一行です．gzのままいけちゃいます．

xr_ensemble = xr.open_dataset('ETOPO1_Ice_g_gmt4.grd.gz')

データセットの編集

系列名を変更:geoviewsで複数の地図をリンクする際に軸名が同じである必要があるため．
描画範囲をスライス：全部は重すぎます．

xr_ensemble = xr_ensemble.rename({'x':'Longitude', 'y':'Latitude'})
xr_ensemble = xr_ensemble.sel(Longitude=slice(95,99), Latitude=slice(14,18))

xarrayをgeoviews.Datasetに変更

この形式以外でも描けますが，geoviews.Datasetは高速なので，データが多い場合は有効です．

kdims = ['Longitude', 'Latitude']
vdims = ['z']
xr_dataset = gv.Dataset(xr_ensemble, kdims=kdims, vdims=vdims)

レンダリング

今回はESRI mapとリンクさせてます．

%%output size=100 filename='out' fig='html'
%opts Image {+framewise} 
xr_dataset.to(gv.Image, ['Longitude', 'Latitude']).redim.range(z=(-50, 10)).options(cmap='jet',colorbar=True,width=450, height=400) \
+ gts.ESRI.options(width=400, height=400) * gts.StamenLabels.options(level='annotation')