またまた地図ネタです．
先日ネタにでましたミャンマーのある河川の河口部について，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')
ソースコードは以下です．
github.com
ETOPOメモ
備忘録です．
ETOPO 5/2/1について
データセットの解像度を示す.ETOPO1が最新.
- ETOPO5 : 5分 緯度5分=~9.2km
- ETOPO2 : 2分 緯度2分=~3.7km
- ETOPO1 : 1分 緯度1分=~1.8km
http://ofgs.aori.u-tokyo.ac.jp/~okino/gmtscripts/datainfo.html
Ice Surface or Bedrock
https://www.ngdc.noaa.gov/mgg/global/etopo1sources.html
南極とグリーンランドの標高値を氷上か基岩かどちらを使うかの差だけです．基本は，絵的には前者を使ったほうが良いですね．
