Frost days with planetary

frost_days_planetary.py

# %% [markdown]
# https://planetarycomputer.microsoft.com/dataset/cil-gdpcir-cc-by#overview

# %%
# required to locate and authenticate with the stac collection
import planetary_computer
import pystac_client

# required to load a zarr array using xarray
import xarray as xr

# optional imports used in this notebook
import pandas as pd
from dask.diagnostics import ProgressBar
from tqdm.auto import tqdm

# %%
catalog = pystac_client.Client.open(
    "https://planetarycomputer.microsoft.com/api/stac/v1",
    modifier=planetary_computer.sign_inplace,
)
collection_cc0 = catalog.get_collection("cil-gdpcir-cc-by")
items = catalog.search(
    collections=["cil-gdpcir-cc-by"],
    query={
        # "cmip6:source_id": {"eq": "FGOALS-g3"},
        "cmip6:experiment_id": {"in": ["historical", "ssp370"]},
    },
)


# %%
# collection_cc0.summaries.to_dict()

# %%
ensemble = items.item_collection()
len(ensemble)


# %%
import collections

collections.Counter(x.collection_id for x in ensemble)

# %%
# select this variable ID for all models in a collection
variable_id = "tasmin"

datasets_by_model = []

for item in tqdm(ensemble):
    asset = item.assets[variable_id]
    datasets_by_model.append(
        xr.open_dataset(asset.href, **asset.extra_fields["xarray:open_kwargs"])
    )

all_datasets = xr.concat(
    datasets_by_model,
    dim=pd.Index([ds.attrs["source_id"] for ds in datasets_by_model], name="model"),
    combine_attrs="drop_conflicts",
)

# %%
all_datasets

# %%
# let's select a subset of the data for the first five days of 2020 over Japan.
# Thanks to https://gist.github.com/graydon/11198540 for the bounding box!
subset = all_datasets.tasmin.sel(
    lon=slice(129.408463169, 145.543137242),
    lat=slice(31.0295791692, 45.5514834662),
    time=slice("2020-01-01", "2022-01-01"),
)

# %%
with ProgressBar():
    subset = subset.compute()

# %%
subset

# %%
subset.mean(dim = 'model').mean(dim = 'time').plot()

# %%
import xclim

frost_days = xclim.indicators.atmos.frost_days(tasmin=subset, thresh='0 degC')
frost_days

# %%
frost_days.mean(dim = 'model').isel(time = 1).plot()