NOAAの海洋熱波カテゴリーをGeoViewsでインタラクティブにプロットする

はじめに

NOAAは水平分解約5kmの海洋熱波カテゴリーを計算したデータを公開している。海洋熱波の強さを1から4の値にした指標である(0は発生なし、-1は海氷)。サイトで見られる図は全球で、例えば日本周辺の細かい様子は見にくい。そこでGeoViewsというPythonのライブラリーでズームなどがインタラクティブにできるプロットにしてみる。

GeoViewsはHoloViewsという可視化を助けてくれるライブラリを、地理的な図も簡単に作れるように拡張したライブラリーである。

コード

作成したJupyter notebookはGithubのgistに置いた。内容は以下の通り。

ライブラリを読み込む。

import geoviews as gv
import geoviews.feature as gf
import xarray as xr
from cartopy import crs
import holoviews as hv
from holoviews.operation.datashader import rasterize
gv.extension('bokeh')

データはnetcdfで、ローカルディスクにダウンロードしている。ここでは2021年のデータを使う。

!ls category/nc/noaa-crw_mhw_v1.0_category_2021*.nc | tail

結果

category/nc/noaa-crw_mhw_v1.0_category_20210806.nc
category/nc/noaa-crw_mhw_v1.0_category_20210807.nc
category/nc/noaa-crw_mhw_v1.0_category_20210808.nc
category/nc/noaa-crw_mhw_v1.0_category_20210809.nc
category/nc/noaa-crw_mhw_v1.0_category_20210810.nc
category/nc/noaa-crw_mhw_v1.0_category_20210811.nc
category/nc/noaa-crw_mhw_v1.0_category_20210812.nc
category/nc/noaa-crw_mhw_v1.0_category_20210813.nc
category/nc/noaa-crw_mhw_v1.0_category_20210814.nc
category/nc/noaa-crw_mhw_v1.0_category_20210815.nc

xarrayを使って読み込む。

ds=xr.open_mfdataset("category/nc/noaa-crw_mhw_v1.0_category_2021*.nc",concat_dim="time",parallel=True)

print(ds)

結果

Dimensions: (time: 226, lat: 3600, lon: 7200)
Coordinates:
* time (time) datetime64[ns] 2021-01-01T12:00:00 ... 2021-08-...
Dimensions without coordinates: lat, lon
Data variables:
latitude (time, lat) float32 dask.array
longitude (time, lon) float32 dask.array
crs (time) int16 -32767 -32767 -32767 ... -32767 -32767
heatwave_category (time, lat, lon) float32 dask.array
mask (time, lat, lon) float32 dask.array
Attributes: (12/58)
Conventions: CF-1.6, ACDD-1.3
ncei_template_version: NCEI_NetCDF_Grid_Template_v2.0
title: NOAA Coral Reef Watch Daily Global 5km Satell...
summary: This is a product of NOAA Coral Reef Watch Da...
references: https://coralreefwatch.noaa.gov/ and https://...
institution: NOAA/NESDIS/STAR Coral Reef Watch program
... ...
publisher_url: https://coralreefwatch.noaa.gov
publisher_email: [e-mail]
contributor_name: NOAA Coral Reef Watch program
contributor_role: Collecting source data and deriving products;...
processing_level: Derived from L4 satellite sea surface tempera...
cdm_data_type: Grid

変な形のデータになっている(経度latitudeが時間の関数の関数になっていたりする)ので、あらためてxarrayのデータセットds2を定義する。heatwave_categoryがプロットしたい変数である。

lon=ds.longitude.isel(time=0)
lat=ds.latitude.isel(time=0)
ds2 = xr.Dataset(
    data_vars=dict(
        heatwave_category=(["time","latitude", "longitude"],ds["heatwave_category"].data)
    ),
    coords=dict(
        time=ds.time.data,
        longitude= lon.data,
        latitude= lat.data,
    )
)

GeoViewsで使えるデータ型に変換する。

kdims = ['latitude','longitude',"time"]
vdims = ['heatwave_category']
hv_dataset = gv.Dataset(ds2, kdims=kdims, vdims=vdims)

GeoViewsで経度緯度座標系のデータをquadmeshのプロットにする。使われていないtimeはセレクターになる。

qmesh=hv_dataset.to(gv.QuadMesh,["longitude","latitude"],crs=crs.PlateCarree()).opts(colorbar=True)

このままだと緯度7200 × 経度3600 x 日数の膨大なデータをインタラクティブに扱うことになる。そこで図を適宜ラスタライズしてくれるライブラリーDatashaderを使ってrasterizeする¹。座標系は東経180°を中心とするように座標変換する。x軸やy軸に緯度経度を表示したほうが便利だが、GeoViewsにはまだバグがあるようなので(central_longitude doesn't update lon tick labels · Issue #329 )、ここではxaxis=None, yaxis=Noneとしておく。

colors = ['gray', 'lightblue', 'yellow', 'orange', 'red', 'darkred']
color_keys = {-1:'gray', 1:'lightblue', 2:'yellow', 3:'orange', 4:'red',5:"darkred"}
raster_qmesh=rasterize(qmesh).opts(clim=(-1, 5),cmap=colors,color_levels=[-1,0,1,2,3,4,5],colorbar=True,projection=crs.PlateCarree(central_longitude=180),global_extent=True,xaxis=None,yaxis=None)

沿岸線と陸地のプロットを用意しておく。

coastline=gf.coastline().opts(scale="10m")
land=gf.land().opts(scale="10m")

ラスター、陸地、沿岸線のプロットを重ねて完成。

(raster_qmesh*land*coastline).opts(width=1000, height=600) 

結果

結果の図は動画のようになる。マウス操作で拡大縮小したり、スライダーで日にちを選んだりできる。ものすごくサクサク動くというわけでもないが、大量のデータを削ることなくインタラクティブにプロットできるという点では十分¹。

---

1. Holoviews(Geoviews)で大量のデータを扱う場合については、Working with large data using datashader — HoloViews documentationを参照。これを参考にDatashaderを使った。文章中のgv.operation.projectを使えばさらに動作が速くなるかもしれないが、うまい使い方を見つけられていない。GPUを使う方法についても書いてあるが、GPUにのせるにはデータが大量すぎるようで、これも今の所うまくいっていない。 ↩