Skip to content

Commit

Permalink
add methane
Browse files Browse the repository at this point in the history
  • Loading branch information
@malmans2
malmans2 committed Mar 26, 2024
1 parent ca23bfd commit 12bb0a4
Showing 1 changed file with 54 additions and 26 deletions.
80 changes: 54 additions & 26 deletions notebooks/wp5/xch4_xco2_satellite_lev3.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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",
Expand All @@ -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)"
]
},
{
Expand All @@ -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",
Expand All @@ -94,79 +106,95 @@
"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",
" da = da.groupby(\"year\").map(diagnostics.seasonal_weighted_mean)\n",
" 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",
Expand Down Expand Up @@ -301,7 +329,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.10.11"
"version": "3.11.7"
}
},
"nbformat": 4,
Expand Down

0 comments on commit 12bb0a4

Please sign in to comment.