{"notebook_content":"\n\n\n\n\nNotebook\n\n\n\n\n\n\n\n\n\n\n\n\n
\n
\n
\n
\n
\n
\n
\n\"No\n

Trends and climatology of NO2\u00b6

\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n

Daily, monthly and annual mean concentrations of NO2 are calculated from one ACTRIS site with a long time series to calculate trends and climatology

\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n

Import libraries\u00b6

\n
\n
\n
\n
\n
\n
\n
\n
\n
In\u00a0[1]:
\n
\n
\n
import json\nimport xarray as xr\n\nimport matplotlib.pyplot as plt\nimport pandas as pd\nimport numpy as np\n\nimport pymannkendall as mk # for Mann-Kendall test\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n

Load data\u00b6

In this example the data is selected by searching in the ACTRIS Data Portal and extracting the URLs to their netCDF files in Thredds. The example shows NO2 from Jungraujoch where the time series for the period 2006-2022 is availible from two different files that needs to be combined.

\n
\n
\n
\n
\n
\n
\n
\n
\n
In\u00a0[2]:
\n
\n
\n
# NO2 at Jungfraujoch from 2006-2020\nURL_first=['https://thredds.nilu.no/thredds/dodsC/ebas_doi/VA/U9/WW/VAU9-WW7S.nc']\n# NO2 at Jungfraujoch from 2021-2022\nURL_last=['https://thredds.nilu.no/thredds/dodsC/ebas_doi/DQ/83/EQ/DQ83-EQBE.nc']\n#combine the two URLs\nfile_urls = URL_first + URL_last\nprint(file_urls)\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
['https://thredds.nilu.no/thredds/dodsC/ebas_doi/VA/U9/WW/VAU9-WW7S.nc', 'https://thredds.nilu.no/thredds/dodsC/ebas_doi/DQ/83/EQ/DQ83-EQBE.nc']\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n

Depending on the type of data, the documentation for each file can be found here: https://data.actris.eu/docs

\n

Documentation on the In-Situ file formats (NetCDF) is available here: https://ebas.pages.nilu.no/ebas-io/fileformat_netcdf/index.html

\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n

Combine the two datasets and calculate daily,monthly and annual mean concentrations\u00b6

\n
\n
\n
\n
\n
\n
\n
\n
\n
In\u00a0[3]:
\n
\n
\n
# define the variable of interest and the year of interest\nvariable_interested_standard_name = 'mass_concentration_of_nitrogen_dioxide_expressed_as_nitrogen_in_air'\n\n# define files with the data of interest\ndaily_means_all = []\nmonthly_means_all = []\nannual_means_all = []\n\n\nfor i in range(len(file_urls)):\n    # Open the netCDF file\n    data = xr.open_dataset(file_urls[i])\n    print(i)\n\n    # check if the dataset contains the variable of interest and print out start and end date of the dataset\n    for var in data.data_vars:\n        if 'standard_name' in data[var].attrs and data[var].attrs['standard_name'] == variable_interested_standard_name:\n            print(' File %d, %s, variable: %s %s Period: %s - %s' % (i, data.attrs['ebas_station_code'], \"{:<35s}\".format(var), \"{:<40s}\".format('('+data.attrs['ebas_station_name']+')'), data.time.min().values, data.time.max().values))\n\n            # get the variable values for all the years\n            variable_values = data[var]\n\n            # Calculate daily and monthly means\n            daily_means = variable_values.resample(time='1D').mean('time')               \n            monthly_means = variable_values.resample(time='ME').mean('time')\n            annual_means = variable_values.resample(time='YE').mean('time')    \n            #add the means to the file that collects all the sites\n            daily_means_all.append(daily_means)\n            monthly_means_all.append(monthly_means)\n            annual_means_all.append(annual_means)            \n            break              \n\n    else:\n        print(f\"File {i} does not contain data \")\ndata.close() \n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
0\n File 0, CH0001G, variable: v_nitrogen_dioxide                  (Jungfraujoch)                           Period: 2006-01-01T00:30:00.000000000 - 2020-12-31T23:30:00.000000000\n1\n File 1, CH0001G, variable: v_nitrogen_dioxide_ug_N_per_m3      (Jungfraujoch)                           Period: 2021-01-01T00:30:00.000000000 - 2022-12-31T23:30:00.000000000\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n

Extract the daily, monthly and annual means to separate files\u00b6

\n
\n
\n
\n
\n
\n
\n
\n
\n
In\u00a0[9]:
\n
\n
\n
def concatenate_time_series(data_arrays, is_annual=False):\n    times = []\n    values = []\n    for data_array in data_arrays:\n        if is_annual:\n            # Set the time to January 1st\n            times.extend(data_array.time.values + np.timedelta64(1, 'D') - np.timedelta64(365, 'D'))\n            values.extend(data_array.values)\n        else:\n            times.extend(data_array.time.values)\n            values.extend(data_array.values)\n    return np.array(times), np.array(values)\n\n# Concatenate daily, monthly, and annual means\ndaily_times, daily_values = concatenate_time_series(daily_means_all)\nmonthly_times, monthly_values = concatenate_time_series(monthly_means_all)\nannual_times, annual_values = concatenate_time_series(annual_means_all, is_annual=True)\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n

Plot the timeserie\u00b6

\n
\n
\n
\n
\n
\n
\n
\n
\n
In\u00a0[10]:
\n
\n
\n
fig, ax1 = plt.subplots(figsize=(10, 6))\n\n# Define colors for each type of data\ndaily_color = 'cornflowerblue'\nmonthly_color = 'cornflowerblue'\nannual_color = 'seagreen'\n\n# Plot monthly and annual means on the primary y-axis\nax1.plot(monthly_times, monthly_values, label='monthly mean', linewidth=3, color=monthly_color)\nax1.plot(annual_times, annual_values, label='annual mean', linewidth=4, color=annual_color)\nax1.set_ylabel('Monthly/Annual NO$_2$ conc. (\u00b5gN/m$^3$)')\nax1.set_ylim(0, 0.3)\n\n# Create a secondary y-axis, and lot daily means on the the secondary y-axis\nax2 = ax1.twinx() \nax2.scatter(daily_times, daily_values, label='daily mean', alpha=0.5, s=10, color=daily_color, edgecolors='black')\nax2.set_ylabel('Daily NO$_2$ conc. (\u00b5gN/m$^3$)')\nax2.set_ylim(0, 3)\nax2.set_xlim(daily_times[0], daily_times[-1])\n\n# Combine legends from both axes\nlines_1, labels_1 = ax1.get_legend_handles_labels()\nlines_2, labels_2 = ax2.get_legend_handles_labels()\nax2.legend(lines_1 + lines_2, labels_1 + labels_2, loc='upper right')\n\n\nplt.title('NO$_2$ at Jungfraujoch, 2016-2022 (note different y-axis scales)')\nplt.show()\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n\"No\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n

Make a function that calculates the Mann Kendall trend and Sen-Theil slope:

\n
\n
\n
\n
\n
\n
\n
\n
\n
In\u00a0[11]:
\n
\n
\n
def calculate_trend(data):\n    # Convert data to a pandas Series\n    series = pd.Series(data)\n    \n    # Calculate the Mann-Kendall test statistic and p-value\n    trend_result = mk.original_test(series)\n    p_value = trend_result.p\n    \n    # Calculate the Sen-Theil slope and intercept\n    slope = trend_result.slope\n    intercept = trend_result.intercept \n    \n    return slope, intercept, p_value\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n

Calculate monthly and annual mean trends and plot

\n
\n
\n
\n
\n
\n
\n
\n
\n
In\u00a0[12]:
\n
\n
\n
# Create a DataFrame from annual_times and annual_values\ndf_annual = pd.DataFrame({\n    'time': annual_times,\n    'value': annual_values\n})\n# Create a DataFrame from monthly_times and monthly_values\ndf_monthly = pd.DataFrame({\n    'time': monthly_times,\n    'value': monthly_values\n})\n\n# Calculate the trend for the annual data\ndata_annual = df_annual.value\nslope_annual, intercept_annual, p_value_annual = calculate_trend(data_annual)\n\nyears = df_annual.time\nsen_annual = [slope_annual * x + intercept_annual for x in range(len(years))]\n\n# Calculate the trend for the monthly data\ndata_month = df_monthly.value\nslope_month, intercept_month, p_value_month = calculate_trend(data_month)\n\nmonths = df_monthly.time\nsen_month = [slope_month * x + intercept_month for x in range(len(months))]\n\n# Create subplots\nfig, axs = plt.subplots(1, 2, figsize=(12, 5))\n\n# Annual data subplot\naxs[0].plot(years, data_annual, color='blue')\n# Add the Theil Sen slope and calculate the annual relative trend. \n# Depending on the p-value, the trend line is solid or dashed\nif 0 < p_value_annual < 0.05: \n    label_annual = f'Sen-Theil Slope: {round(100 * slope_annual / intercept_annual, 1)} %/y' \n    axs[0].plot(years, sen_annual, 'darkred', label=label_annual)\nelse:\n    label_annual = f'Sen-Theil Slope: {round(100 * slope_annual / intercept_annual, 1)} %/y'\n    axs[0].plot(years, sen_annual, 'darkred', linestyle='--', label=label_annual)\n\naxs[0].set_yticks(np.arange(0, 0.2, step=0.05))\naxs[0].set_ylabel(\"\u00b5gN/ m$^{3}$\")\naxs[0].set_title('NO$_2$ annual mean trends at Jungfraujoch')\naxs[0].legend()\n\n# Monthly data subplot\naxs[1].plot(months, data_month, color='blue')\n# Add the Theil Sen slope, multiply the monthly Sen slope with 12 to get the annual %trend. \nif 0 < p_value_month < 0.05:\n    label_month = f'Sen-Theil Slope: {round(12 * 100 * slope_month / intercept_month, 1)} %/y' \n    axs[1].plot(months, sen_month, 'darkred', label=label_month)\nelse:\n    label_month = f'Sen-Theil Slope: {round(12 * 100 * slope_month / intercept_month, 1)} %/y'\n    axs[1].plot(months, sen_month, 'darkred', linestyle='--', label=label_month)\n\n\naxs[1].set_yticks(np.arange(0, 0.4, step=0.1))\naxs[1].set_ylabel(\"\u00b5gN/ m$^{3}$\")\naxs[1].set_title('NO$_2$ monthly mean trends at Jungfraujoch')\naxs[1].legend()\n\nplt.tight_layout()\nplt.show()\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n\"No\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n

Plot monthly mean of 2022 and compare with earlier period (climatology)\u00b6

\n
\n
\n
\n
\n
\n
\n
\n
\n
In\u00a0[13]:
\n
\n
\n
# Extract the month from the time column\ndf_monthly['month'] = df_monthly['time'].dt.month\n\n# Filter the DataFrame to create two separate DataFrames,\n# one with 2022 data and one with other years\ndf_2022 = df_monthly[df_monthly['time'].dt.year == 2022].drop(columns=['time'])\ndf_other_years = df_monthly[df_monthly['time'].dt.year != 2022].drop(columns=['time'])\n\n# Calculate monthly statistics for other years and extract relevant statistics\nmonthly_stats = df_other_years.groupby('month')['value'].describe()\nmax_values = monthly_stats['max']\nmin_values = monthly_stats['min']\nq75_values = monthly_stats['75%']\nq25_values = monthly_stats['25%']\nmedian_values = monthly_stats['50%']\nmonths = list(range(1, 13))  # List of integers from 1 to 12\n\n# Plot the statistics\nplt.figure(figsize=(9, 5))\n\nplt.plot(df_2022['month'], df_2022['value'], label='2022', marker='o', color='red')\nplt.plot(months, median_values, label='Median 2016-2021', marker='o')\nplt.plot(months, q75_values, color='black', alpha=0.1)\nplt.plot(months, q25_values, color='black', alpha=0.1)\nplt.plot(months, max_values, color='black', alpha=0.2)\nplt.plot(months, min_values, color='black', alpha=0.2)\n\n# Add shaded area between Max and 75%\nplt.fill_between(months, q25_values, q75_values, color='gray', alpha=0.6, label='25% - 75% percentiles, 2016-2021')\nplt.fill_between(months, q75_values, max_values, color='gray', alpha=0.3,label='75% (25%) percentiles - max (min), 2016-2021')\n\n# Add shaded area between Max and 75%\nplt.fill_between(months, q25_values, min_values, color='gray', alpha=0.3)\n\nplt.xlim(1, 12)\nplt.yticks(np.arange(0, 0.4, 0.1))\nplt.xticks(months, ['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec'])\nplt.ylabel('NO$_2$ conc. (\u00b5gN/m$^3$)')\nplt.title('Climatology of monthly mean concentrations of NO$_2$ at Jungfraujoch, 2016-2022')\nplt.legend()\nplt.grid(True, linestyle='--', alpha=0.5)\nplt.show()\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n
\n\"No\n
\n
\n
\n
\n
\n
\n\n\n"}