趣味で計算流砂水理

Computational Sediment Hydraulics for Fun

ミャンマー某河川の河口部:ETOPO1 + geoviews

またまた地図ネタです.

先日ネタにでましたミャンマーのある河川の河口部について,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

南極とグリーンランドの標高値を氷上か基岩かどちらを使うかの差だけです.基本は,絵的には前者を使ったほうが良いですね.