From 12bb0a491f47cf74ca076ce79199abcda90e2c13 Mon Sep 17 00:00:00 2001 From: Mattia Almansi Date: Tue, 26 Mar 2024 12:57:19 +0100 Subject: [PATCH] add methane --- notebooks/wp5/xch4_xco2_satellite_lev3.ipynb | 80 +++++++++++++------- 1 file changed, 54 insertions(+), 26 deletions(-) diff --git a/notebooks/wp5/xch4_xco2_satellite_lev3.ipynb b/notebooks/wp5/xch4_xco2_satellite_lev3.ipynb index 307ab41..bb2bed8 100644 --- a/notebooks/wp5/xch4_xco2_satellite_lev3.ipynb +++ b/notebooks/wp5/xch4_xco2_satellite_lev3.ipynb @@ -25,7 +25,7 @@ "import cartopy.crs as ccrs\n", "import matplotlib.pyplot as plt\n", "import xarray as xr\n", - "from c3s_eqc_automatic_quality_control import diagnostics, download, plot\n", + "from c3s_eqc_automatic_quality_control import diagnostics, download, plot, utils\n", "from xarrayMannKendall import Mann_Kendall_test\n", "\n", "plt.style.use(\"seaborn-v0_8-notebook\")\n", @@ -48,10 +48,22 @@ "source": [ "# Choose variable (xch4 or xco2)\n", "variable = \"xco2\"\n", - "assert variable in (\"xch4\", \"xco2\")\n", + "assert variable in [\n", + " f\"{prefix}{suffix}\"\n", + " for prefix in (\"xch4\", \"xco2\")\n", + " for suffix in (\"\", \"_nobs\", \"_stderr\", \"_stddev\")\n", + "]\n", "\n", "# Minimum value of land fraction used for masking\n", - "min_land_fraction = None # None: Do not apply mask" + "min_land_fraction = 0.5 # None: Do not apply mask\n", + "\n", + "# Choose a time period\n", + "year_start = 2012\n", + "year_stop = 2021\n", + "\n", + "# Define region for analysis\n", + "lon_slice = slice(-180, 180)\n", + "lat_slice = slice(-90, 90)" ] }, { @@ -69,10 +81,10 @@ "outputs": [], "source": [ "request = (\n", - " \"satellite-carbon-dioxide\" if variable == \"xco2\" else \"satellite-methane\",\n", + " \"satellite-carbon-dioxide\" if variable.startswith(\"xco2\") else \"satellite-methane\",\n", " {\n", " \"processing_level\": \"level_3\",\n", - " \"variable\": variable,\n", + " \"variable\": variable.split(\"_\")[0],\n", " \"sensor_and_algorithm\": \"merged_obs4mips\",\n", " \"version\": \"4.4\",\n", " \"format\": \"zip\",\n", @@ -94,31 +106,43 @@ "metadata": {}, "outputs": [], "source": [ - "transform_func_kwargs = {\"min_land_fraction\": min_land_fraction}\n", + "transform_func_kwargs = {\n", + " \"min_land_fraction\": min_land_fraction,\n", + " \"variable\": variable,\n", + " \"year_start\": year_start,\n", + " \"year_stop\": year_stop,\n", + " \"lon_slice\": lon_slice,\n", + " \"lat_slice\": lat_slice,\n", + "}\n", "\n", "\n", - "def get_da(ds, min_land_fraction):\n", - " (varname,) = set(ds.data_vars) & {\"xch4\", \"xco2\"}\n", - " da = ds[varname]\n", + "def get_da(\n", + " ds, min_land_fraction, variable, year_start, year_stop, lon_slice, lat_slice\n", + "):\n", + " da = ds[variable].sel(time=slice(str(year_start), str(year_stop)))\n", + " da = utils.regionalise(da, lon_slice=lon_slice, lat_slice=lat_slice)\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", + " if da.name.endswith(\"_nobs\"):\n", + " return da\n", + "\n", " with xr.set_options(keep_attrs=True):\n", - " if da.name == \"xch4\" and da.attrs[\"units\"] != \"ppb\":\n", + " if da.name.startswith(\"xch4\") and da.attrs[\"units\"] != \"ppb\":\n", " da *= 1.0e9\n", " da.attrs[\"units\"] = \"ppb\"\n", - " elif da.name == \"xco2\" and da.attrs[\"units\"] != \"ppm\":\n", + " elif da.name.startswith(\"xco2\") and da.attrs[\"units\"] != \"ppm\":\n", " da *= 1.0e6\n", " da.attrs[\"units\"] = \"ppm\"\n", " return da\n", "\n", "\n", - "def compute_seasonal_timeseries(ds, min_land_fraction):\n", + "def compute_seasonal_timeseries(ds, **get_da_kwargs):\n", " # Shift years (shift -1 to get D(year-1)J(year)F(year))\n", - " da = get_da(ds, min_land_fraction)\n", + " da = get_da(ds, **get_da_kwargs)\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", @@ -126,47 +150,51 @@ " return convert_units(da).to_dataset()\n", "\n", "\n", - "def compute_statistics(ds, min_land_fraction):\n", - " da = get_da(ds, min_land_fraction)\n", + "def compute_statistics(ds, **get_da_kwargs):\n", + " da = get_da(ds, **get_da_kwargs)\n", " da = diagnostics.spatial_weighted_statistics(da)\n", " return convert_units(da).to_dataset()\n", "\n", "\n", - "def compute_monthly_anomalies(ds, min_land_fraction):\n", - " da = get_da(ds, min_land_fraction)\n", + "def compute_monthly_anomalies(ds, **get_da_kwargs):\n", + " da = get_da(ds, **get_da_kwargs)\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", "\n", "\n", - "def compute_mann_kendall_trend(da, **kwargs):\n", + "def compute_mann_kendall_trend(da, **mann_kendall_kwargs):\n", " coords_name = {\"time\": \"time\", \"x\": \"longitude\", \"y\": \"latitude\"}\n", - " ds_trend = Mann_Kendall_test(da, coords_name=coords_name, **kwargs).compute()\n", + " ds_trend = Mann_Kendall_test(\n", + " da, coords_name=coords_name, **mann_kendall_kwargs\n", + " ).compute()\n", " return ds_trend.rename({k: v for k, v in coords_name.items() if k != \"time\"})\n", "\n", "\n", - "def compute_seasonal_detrended_anomaly(da, **kwargs):\n", - " da_trend = xr.polyval(da[\"time\"], da.polyfit(\"time\", **kwargs).polyfit_coefficients)\n", + "def compute_seasonal_detrended_anomaly(da, **polyfit_kwargs):\n", + " da_trend = xr.polyval(\n", + " da[\"time\"], da.polyfit(\"time\", **polyfit_kwargs).polyfit_coefficients\n", + " )\n", " da_detrended = da - da_trend\n", " return da_detrended.groupby(\"time.year\").map(diagnostics.seasonal_weighted_mean)\n", "\n", "\n", - "def compute_anomaly_trends(ds, min_land_fraction):\n", - " da_anomaly = compute_monthly_anomalies(ds, min_land_fraction)\n", + "def compute_anomaly_trends(ds, **get_da_kwargs):\n", + " da_anomaly = compute_monthly_anomalies(ds, **get_da_kwargs)\n", "\n", " # Mann-Kendall\n", " ds_mann_kendall = compute_mann_kendall_trend(\n", " da_anomaly, alpha=0.05, method=\"theilslopes\"\n", " ).where(da_anomaly.notnull().any(\"time\"))\n", " ds_mann_kendall[\"trend\"].attrs = {\n", - " \"long_name\": f\"Trend of anomalies of {da_anomaly.attrs['long_name']}\",\n", + " \"long_name\": f\"Trend of anomalies of {da_anomaly.attrs.get('long_name', da_anomaly.name)}\",\n", " \"units\": f\"{da_anomaly.attrs['units']}/month\",\n", " }\n", "\n", " # Detrended anomalies\n", " da_detrended = compute_seasonal_detrended_anomaly(da_anomaly, deg=1)\n", " da_detrended.attrs = {\n", - " \"long_name\": f\"Detrended anomalies of {da_anomaly.attrs['long_name']}\",\n", + " \"long_name\": f\"Detrended of anomalies of {da_anomaly.attrs.get('long_name', da_anomaly.name)}\",\n", " \"units\": f\"{da_anomaly.attrs['units']}\",\n", " }\n", "\n", @@ -301,7 +329,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.11" + "version": "3.11.7" } }, "nbformat": 4,