Skip to content

Commit

Permalink
new esacci
Browse files Browse the repository at this point in the history
  • Loading branch information
@malmans2
malmans2 committed Aug 11, 2023
1 parent 753ab9c commit b3b1541
Showing 1 changed file with 145 additions and 105 deletions.
250 changes: 145 additions & 105 deletions notebooks/wp5/satellite_esacci_gmpe.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -27,6 +27,7 @@
"source": [
"import cartopy.crs as ccrs\n",
"import matplotlib.pyplot as plt\n",
"import tqdm\n",
"import xarray as xr\n",
"from c3s_eqc_automatic_quality_control import diagnostics, download, plot, utils\n",
"\n",
Expand All @@ -50,8 +51,8 @@
"outputs": [],
"source": [
"# Time\n",
"start = \"1982-01\"\n",
"stop = \"1991-12\"\n",
"year_start = 1982\n",
"year_stop = 2011\n",
"\n",
"# Regions\n",
"regions = {\n",
Expand All @@ -78,7 +79,7 @@
"source": [
"# Requests\n",
"request_dicts = {\n",
" \"esacci\": {\n",
" \"ESACCI\": {\n",
" \"collection_id\": \"satellite-sea-surface-temperature\",\n",
" \"request\": {\n",
" \"processinglevel\": \"level_4\",\n",
Expand All @@ -89,7 +90,7 @@
" },\n",
" \"chunks\": {\"year\": 1, \"month\": 1},\n",
" },\n",
" \"gmpe\": {\n",
" \"GMPE\": {\n",
" \"collection_id\": \"satellite-sea-surface-temperature-ensemble-product\",\n",
" \"request\": {\n",
" \"format\": \"zip\",\n",
Expand Down Expand Up @@ -130,53 +131,44 @@
" da = ds[\"analysed_sst\"]\n",
" if \"mask\" in ds:\n",
" da = da.where(ds[\"mask\"] == 1)\n",
" with xr.set_options(keep_attrs=True):\n",
" da -= 273.15\n",
" da.attrs[\"units\"] = \"°C\"\n",
" return da\n",
"\n",
"\n",
"def rechunk(obj):\n",
" \"\"\"Use NetCDF chunks.\"\"\"\n",
" chunks = {\"season\": 1, \"latitude\": 1_200, \"longitude\": 2_400}\n",
" chunks = {\"year\": 1, \"season\": 1, \"latitude\": 1_200, \"longitude\": 2_400}\n",
" return obj.chunk(\n",
" **{dim: chunksize for dim, chunksize in chunks.items() if dim in obj.dims}\n",
" )\n",
"\n",
"\n",
"def compute_time_reductions(ds, reductions, groups):\n",
"def compute_time_reductions(ds, func, **kwargs):\n",
" ds = rechunk(ds)\n",
" da = get_masked_sst(ds)\n",
" dataarrays = []\n",
" for group in groups:\n",
" for reduction in reductions:\n",
" name = \"_\".join([group, reduction])\n",
" long_name = f\"{reduction.title()} of {da.attrs['long_name']}\"\n",
" func = getattr(diagnostics, f\"{group}_weighted_{reduction}\")\n",
" da_reduced = rechunk(func(da, weights=False))\n",
" da_reduced.attrs[\"long_name\"] = long_name\n",
" da_reduced.encoding[\"chunksizes\"] = tuple(map(max, da_reduced.chunks))\n",
" dataarrays.append(da_reduced.rename(name))\n",
" return xr.merge(dataarrays)\n",
" da_reduced = rechunk(func(da, **kwargs))\n",
" if \"season\" in da_reduced.dims:\n",
" da_reduced = da_reduced.sel(season=sorted(set(da[\"time\"].dt.season.values)))\n",
" da_reduced.encoding[\"chunksizes\"] = tuple(map(max, da_reduced.chunks))\n",
" return da_reduced.to_dataset()\n",
"\n",
"\n",
"def compute_spatial_weighted_reductions(ds, reductions, lon_slice, lat_slice):\n",
"def compute_spatial_weighted_reductions(ds, reduction, lon_slice, lat_slice):\n",
" ds = rechunk(ds)\n",
" ds = utils.regionalise(ds, lon_slice=lon_slice, lat_slice=lat_slice)\n",
" da = get_masked_sst(ds)\n",
" da = diagnostics.spatial_weighted_mean(da, weights=True).compute()\n",
" grouped = da.groupby(\"time.dayofyear\")\n",
" dataarrays = []\n",
" for reduction in reductions:\n",
" long_name = f\"{reduction.title()} of {da.attrs['long_name']}\"\n",
" func = getattr(grouped, reduction)\n",
" da_reduced = rechunk(func(\"time\", keep_attrs=True))\n",
" da_reduced.attrs[\"long_name\"] = long_name\n",
" dataarrays.append(da_reduced.rename(reduction))\n",
" return xr.merge(dataarrays)"
" func = getattr(grouped, reduction)\n",
" da_reduced = func(\"time\", keep_attrs=True)\n",
" return da_reduced.to_dataset()"
]
},
{
"attachments": {},
"cell_type": "markdown",
"id": "93ad0c07",
"id": "59ddd578-1703-47f3-bc84-6d4d8b877a11",
"metadata": {},
"source": [
"## Download and transform"
Expand All @@ -185,111 +177,159 @@
{
"cell_type": "code",
"execution_count": null,
"id": "2e1742e5-88f7-465b-84c2-6b8e770847e8",
"id": "6c2a4c0a-b185-44cb-9eca-26cde95da101",
"metadata": {},
"outputs": [],
"source": [
"datasets_maps = {}\n",
"datasets_timeseries = []\n",
"\n",
"# Settings\n",
"reductions = (\"mean\", \"std\")\n",
"groups = (\"time\", \"seasonal\")\n",
"season_month_dict = {\n",
" \"DJF\": {\"12\", \"01\", \"02\"},\n",
" \"MAM\": {\"03\", \"04\", \"05\"},\n",
" \"JJA\": {\"06\", \"07\", \"08\"},\n",
" \"SON\": {\"09\", \"10\", \"11\"},\n",
"}\n",
"\n",
"# Initialize variables\n",
"datasets_annual = {}\n",
"datasets_seasonal = {}\n",
"timeseries_ds_list = []\n",
"for product, request_dict in request_dicts.items():\n",
" # Common kwargs\n",
" kwargs = {\n",
" \"collection_id\": request_dict[\"collection_id\"],\n",
" \"requests\": download.update_request_date(\n",
" request_dict[\"request\"], start=start, stop=stop, stringify_dates=True\n",
" ),\n",
" \"transform_chunks\": False,\n",
" \"chunks\": request_dict[\"chunks\"],\n",
" **open_mfdataset_kwargs,\n",
" }\n",
"\n",
" # Time reductions\n",
" ds = download.download_and_transform(\n",
" **kwargs,\n",
" transform_func=compute_time_reductions,\n",
" transform_func_kwargs={\"reductions\": reductions, \"groups\": groups},\n",
" # Annual\n",
" print(f\"{product=} {reductions=}\")\n",
" for year in tqdm.tqdm(range(year_start, year_stop + 1), desc=\"annual\"):\n",
" annual_requests = download.update_request_date(\n",
" request_dict[\"request\"],\n",
" start=f\"{year}-01\",\n",
" stop=f\"{year}-12\",\n",
" stringify_dates=True,\n",
" )\n",
" annual_ds_list = []\n",
" for reduction in reductions:\n",
" func = getattr(diagnostics, f\"annual_weighted_{reduction}\")\n",
" ds = download.download_and_transform(\n",
" **kwargs,\n",
" requests=annual_requests,\n",
" transform_func=compute_time_reductions,\n",
" transform_func_kwargs={\"func\": func, \"weights\": False},\n",
" )\n",
" annual_ds_list.append(rechunk(ds).expand_dims(reduction=[reduction]))\n",
" datasets_annual[product] = xr.concat(annual_ds_list, \"reduction\")\n",
"\n",
" # Seasonal and spatial\n",
" requests = download.update_request_date(\n",
" request_dict[\"request\"],\n",
" start=f\"{year_start}-01\",\n",
" stop=f\"{year_stop}-12\",\n",
" stringify_dates=True,\n",
" )\n",
" datasets_maps[product] = rechunk(ds)\n",
"\n",
" # Spatial weighted reductions\n",
" for region, slices in regions.items():\n",
" ds = download.download_and_transform(\n",
" **kwargs,\n",
" transform_func=compute_spatial_weighted_reductions,\n",
" transform_func_kwargs={\"reductions\": reductions} | slices,\n",
" )\n",
" datasets_timeseries.append(ds.expand_dims(region=[region], product=[product]))\n",
"ds_timeseries = xr.merge(datasets_timeseries)"
]
},
{
"attachments": {},
"cell_type": "markdown",
"id": "460d3bc2",
"metadata": {},
"source": [
"## Plot mean and std maps"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "2d756796-1075-44bf-aadf-d837d369ca89",
"metadata": {},
"outputs": [],
"source": [
"maps_kwargs = {\"projection\": ccrs.Robinson(), \"plot_func\": \"contourf\", \"col_wrap\": 2}\n",
" seasonal_ds_list = []\n",
" for reduction in reductions:\n",
" # Spatial\n",
" print(f\"{product=} {reduction=}\")\n",
" for region, slices in tqdm.tqdm(regions.items(), desc=\"region\"):\n",
" ds = download.download_and_transform(\n",
" **kwargs,\n",
" requests=requests,\n",
" transform_func=compute_spatial_weighted_reductions,\n",
" transform_func_kwargs={\"reduction\": reduction} | slices,\n",
" )\n",
" timeseries_ds_list.append(\n",
" ds.expand_dims(\n",
" product=[product], reduction=[reduction], region=[region]\n",
" )\n",
" )\n",
"\n",
"for product, ds in datasets_maps.items():\n",
" for var, da in ds.data_vars.items():\n",
" plot.projected_map(\n",
" da if \"season\" in da.dims else da.compute(),\n",
" cmap=\"Spectral_r\" if var.endswith(\"mean\") else \"tab20b\",\n",
" col=\"season\" if \"season\" in da.dims else None,\n",
" **maps_kwargs,\n",
" # Seasonal\n",
" func = getattr(diagnostics, f\"seasonal_weighted_{reduction}\")\n",
" tmp_ds_list = []\n",
" for season, months in tqdm.tqdm(season_month_dict.items(), desc=\"season\"):\n",
" season_requests = [\n",
" {\n",
" k: v\n",
" if k != \"month\"\n",
" else sorted(set({v} if isinstance(v, str) else v) & months)\n",
" for k, v in r.items()\n",
" }\n",
" for r in requests\n",
" ]\n",
" ds = download.download_and_transform(\n",
" **kwargs,\n",
" requests=season_requests,\n",
" transform_func=compute_time_reductions,\n",
" transform_func_kwargs={\"func\": func, \"weights\": False},\n",
" )\n",
" tmp_ds_list.append(rechunk(ds))\n",
" seasonal_ds_list.append(\n",
" xr.concat(tmp_ds_list, \"season\").expand_dims(reduction=[reduction])\n",
" )\n",
" plt.suptitle(f\"{product.upper()} ({start}, {stop})\")\n",
" plt.show()"
" datasets_seasonal[product] = xr.concat(seasonal_ds_list, \"reduction\")\n",
"\n",
"ds_timeseries = xr.merge(timeseries_ds_list)\n",
"del tmp_ds_list, annual_ds_list, seasonal_ds_list, timeseries_ds_list"
]
},
{
"attachments": {},
"cell_type": "markdown",
"id": "48d79eb9",
"id": "ed77ac70-113b-4a32-953f-05f8bd1de16c",
"metadata": {},
"source": [
"## Plot bias"
"## Plot low resolution seasonal maps"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "503829f6-0758-4026-ae66-850b66d5143c",
"id": "7e91c57c-a05d-4f21-9184-47ee2085c151",
"metadata": {},
"outputs": [],
"source": [
"for var in (\"time_mean\", \"seasonal_mean\"):\n",
" da = datasets_maps[list(datasets_maps)[0]][var]\n",
"datasets = []\n",
"for product, ds in datasets_seasonal.items():\n",
" if product == \"ESACCI\":\n",
" ds = ds.coarsen(latitude=5, longitude=5).mean()\n",
" ds[\"latitude\"] = ds[\"latitude\"].round(3)\n",
" ds[\"longitude\"] = ds[\"longitude\"].round(3)\n",
" datasets.append(ds.expand_dims(product=[product]))\n",
"ds = xr.concat(datasets, \"product\")\n",
"\n",
"for reduction, da in ds[\"analysed_sst\"].groupby(\"reduction\"):\n",
" title = f\"{reduction.title()} ({year_start}-{year_stop})\"\n",
" kwargs = {\"projection\": ccrs.Robinson(), \"plot_func\": \"contourf\"}\n",
"\n",
" # Reduction\n",
" plot.projected_map(\n",
" da,\n",
" cmap=\"Spectral_r\" if reduction == \"mean\" else \"tab20b\",\n",
" center=False,\n",
" col=\"product\",\n",
" row=\"season\",\n",
" **kwargs,\n",
" )\n",
" plt.suptitle(title)\n",
" plt.show()\n",
"\n",
" # Bias\n",
" with xr.set_options(keep_attrs=True):\n",
" da = da - datasets_maps[list(datasets_maps)[1]][var].interp_like(da)\n",
" da.attrs[\"long_name\"] = da.attrs[\"long_name\"].replace(\"Mean\", \"Mean bias\")\n",
" bias = da.diff(\"product\").squeeze()\n",
" plot.projected_map(\n",
" da if \"season\" in da.dims else da.compute(),\n",
" cmap=\"PRGn\",\n",
" col=\"season\" if \"season\" in da.dims else None,\n",
" **maps_kwargs,\n",
" bias, cmap=\"PRGn\", center=True, col=\"season\", col_wrap=2, **kwargs\n",
" )\n",
" plt.suptitle(f\"{' - '.join(list(datasets_maps)).upper()} ({start}, {stop})\")\n",
" plt.suptitle(\"Bias of \" + title + \"\\n\" + \" - \".join(da[\"product\"][::-1].values))\n",
" plt.show()"
]
},
{
"attachments": {},
"cell_type": "markdown",
"id": "b1a309a7",
"id": "13fac4cb-afc6-411a-9dc8-497423d8e900",
"metadata": {},
"source": [
"## Plot timeseries"
Expand All @@ -298,25 +338,25 @@
{
"cell_type": "code",
"execution_count": null,
"id": "626e80e5-f301-4a11-9ff4-52dbe5c71815",
"id": "5fda07b8-221a-47a2-a8bd-c7ad133bab17",
"metadata": {},
"outputs": [],
"source": [
"for region, ds_region in ds_timeseries.groupby(\"region\"):\n",
"for region, da_region in ds_timeseries[\"analysed_sst\"].groupby(\"region\"):\n",
" fig, ax = plt.subplots()\n",
" ax.set_prop_cycle(color=[\"green\", \"blue\"])\n",
" for product, ds_product in ds_region.groupby(\"product\"):\n",
" ds_product[\"mean\"].plot(\n",
" hue=\"product\", ax=ax, label=product.upper(), add_legend=False\n",
" )\n",
" for product, da_product in da_region.groupby(\"product\"):\n",
" mean = da_product.sel(reduction=\"mean\")\n",
" std = da_product.sel(reduction=\"std\")\n",
" mean.plot(hue=\"product\", ax=ax, label=product, add_legend=False)\n",
" ax.fill_between(\n",
" ds_product[\"dayofyear\"],\n",
" ds_product[\"mean\"] - ds_product[\"std\"],\n",
" ds_product[\"mean\"] + ds_product[\"std\"],\n",
" da_product[\"dayofyear\"],\n",
" mean - std,\n",
" mean + std,\n",
" alpha=0.5,\n",
" label=f\"{product.upper()} ± std\",\n",
" )\n",
" ax.set_title(f\"{region.title()} ({start}, {stop})\")\n",
" ax.set_title(f\"{region.title()} ({year_start}-{year_stop})\")\n",
" ax.legend()\n",
" ax.grid()\n",
" plt.show()"
Expand Down

0 comments on commit b3b1541

Please sign in to comment.