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

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

[GeoViews](https://geoviews.org/)は[HoloViews](https://holoviews.org/)という可視化を助けてくれるライブラリを、地理的な図も簡単に作れるように拡張したライブラリーである。

# コード
作成したJupyter notebookはGithubの[gist](https://gist.github.com/tmiyama/08ada0d4b287750117e37a1f161e735f)に置いた。内容は以下の通り。

ライブラリを読み込む。

```python
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年のデータを使う。

```python
!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](https://xarray.pydata.org/en/stable/)を使って読み込む。

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


```python
print(ds)
```

結果
><xarray.Dataset>
    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<chunksize=(1, 3600), meta=np.ndarray>
        longitude          (time, lon) float32 dask.array<chunksize=(1, 7200), meta=np.ndarray>
        crs                (time) int16 -32767 -32767 -32767 ... -32767 -32767
        heatwave_category  (time, lat, lon) float32 dask.array<chunksize=(1, 3600, 7200), meta=np.ndarray>
        mask               (time, lat, lon) float32 dask.array<chunksize=(1, 3600, 7200), meta=np.ndarray>
    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:            coralreefwatch@noaa.gov
        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がプロットしたい変数である。

```python
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で使えるデータ型に変換する。

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

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

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

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

```python
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)
```

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

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

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

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

# 結果

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

<iframe width="560" height="315" src="https://www.youtube.com/embed/BKB55DcKCQ4?rel=0" title="YouTube video player" frameborder="0" allow="accelerometer; autoplay; clipboard-write; encrypted-media; gyroscope; picture-in-picture" allowfullscreen></iframe>

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