add masking

notebooks/wp5/wip_xch4_xco2.ipynb → notebooks/wp5/xch4_xco2_satellite_lev3.ipynb

@@ -48,7 +48,10 @@
    "source": [
     "# Choose variable (xch4 or xco2)\n",
     "variable = \"xco2\"\n",
-    "assert variable in (\"xch4\", \"xco2\")"
+    "assert variable in (\"xch4\", \"xco2\")\n",
+    "\n",
+    "# Minimum value of land fraction used for masking\n",
+    "min_land_fraction = None  # None: Do not apply mask"
    ]
   },
   {
@@ -91,9 +94,14 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "def get_da(ds):\n",
+    "transform_func_kwargs = {\"min_land_fraction\": min_land_fraction}\n",
+    "\n",
+    "def get_da(ds, min_land_fraction):\n",
     "    (varname,) = set(ds.data_vars) & {\"xch4\", \"xco2\"}\n",
-    "    return ds[varname]\n",
+    "    da = ds[varname]\n",
+    "    if min_land_fraction is not None:\n",
+    "        return da.where(ds[\"land_fraction\"] >= min_land_fraction)\n",
+    "    return da\n",
     "\n",
     "\n",
     "def convert_units(da):\n",
@@ -107,24 +115,24 @@
     "    return da\n",
     "\n",
     "\n",
-    "def compute_seasonal_timeseries(ds):\n",
+    "def compute_seasonal_timeseries(ds, min_land_fraction):\n",
     "    # Shift years (shift -1 to get D(year-1)J(year)F(year))\n",
-    "    da = get_da(ds)\n",
+    "    da = get_da(ds, min_land_fraction)\n",
     "    da = da.assign_coords(year=ds[\"time\"].dt.year.shift(time=-1).astype(int))\n",
     "    # Get rid of 1st JF and last D, so it become [MAM, JJA, SON, DJF, ..., SON]\n",
     "    da = da.isel(time=slice(2, -1))\n",
     "    da = da.groupby(\"year\").map(diagnostics.seasonal_weighted_mean)\n",
     "    return convert_units(da).to_dataset()\n",
     "\n",
     "\n",
-    "def compute_statistics(ds):\n",
-    "    da = get_da(ds)\n",
+    "def compute_statistics(ds, min_land_fraction):\n",
+    "    da = get_da(ds, min_land_fraction)\n",
     "    da = diagnostics.spatial_weighted_statistics(da)\n",
     "    return convert_units(da).to_dataset()\n",
     "\n",
     "\n",
-    "def compute_monthly_anomalies(ds):\n",
-    "    da = get_da(ds)\n",
+    "def compute_monthly_anomalies(ds, min_land_fraction):\n",
+    "    da = get_da(ds, min_land_fraction)\n",
     "    with xr.set_options(keep_attrs=True):\n",
     "        da = da.groupby(\"time.month\") - da.groupby(\"time.month\").mean()\n",
     "    return convert_units(da)\n",
@@ -142,8 +150,8 @@
     "    return da_detrended.groupby(\"time.year\").map(diagnostics.seasonal_weighted_mean)\n",
     "\n",
     "\n",
-    "def compute_trends(ds):\n",
-    "    da_anomaly = compute_monthly_anomalies(ds)\n",
+    "def compute_trends(ds, min_land_fraction):\n",
+    "    da_anomaly = compute_monthly_anomalies(ds, min_land_fraction)\n",
     "\n",
     "    # Mann-Kendall\n",
     "    ds_mann_kendall = compute_mann_kendall_trend(\n",
@@ -180,7 +188,9 @@
    "outputs": [],
    "source": [
     "ds_seasonal = download.download_and_transform(\n",
-    "    *request, transform_func=compute_seasonal_timeseries\n",
+    "    *request,\n",
+    "    transform_func=compute_seasonal_timeseries,\n",
+    "    transform_func_kwargs=transform_func_kwargs,\n",
     ")\n",
     "_ = plot.projected_map(\n",
     "    ds_seasonal[variable],\n",
@@ -205,7 +215,11 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "ds_stats = download.download_and_transform(*request, transform_func=compute_statistics)\n",
+    "ds_stats = download.download_and_transform(\n",
+    "    *request,\n",
+    "    transform_func=compute_statistics,\n",
+    "    transform_func_kwargs=transform_func_kwargs,\n",
+    ")\n",
     "fig, ax = plt.subplots(1, 1)\n",
     "ds_stats[variable].drop_sel(diagnostic=\"std\").plot(hue=\"diagnostic\", ax=ax)\n",
     "mean = ds_stats[variable].sel(diagnostic=\"mean\")\n",
@@ -228,7 +242,9 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "ds_trend = download.download_and_transform(*request, transform_func=compute_trends)\n",
+    "ds_trend = download.download_and_transform(\n",
+    "    *request, transform_func=compute_trends, transform_func_kwargs=transform_func_kwargs\n",
+    ")\n",
     "\n",
     "plot.projected_map(ds_trend[\"trend\"], robust=True, projection=ccrs.Robinson())\n",
     "plot.projected_map(\n",