\n $\"No$ \n

Daily, monthly and annual mean concentrations of NO2\u00b6

Example script to plot the timeseries of daily and monthly mean of NO2 mass concentration data for 2022 at ACTRIS stations. Furher, the annual means are plotted on a European map.

Import libraries\u00b6

In\u00a0[2]:

import xarray as xr\nimport matplotlib.pyplot as plt\nimport matplotlib as mpl\nimport pandas as pd\nimport numpy as np\nimport json\nimport cartopy.crs as ccrs\nimport cartopy.feature as cfeature\nimport cartopy.io.img_tiles as cimgt\nfrom matplotlib.colors import BoundaryNorm\nfrom matplotlib import colormaps as cm\nimport os\n

Load data\u00b6

In this example the data is selected by searching in the ACTRIS Data Portal. A basket (.json) file from choosing variable Nitrogen dioxide mass concentration and framework ACTRIS is stored in static. A more thorugh example of this process is shown in file \"Search and access/Search with the ACTRIS Data Portal.ipynb\".

In\u00a0[3]:

os.chdir(os.getcwd())\n#named the basket file downloaded from ACTRIS Data Portal \"basket_no2_2022.json\"\nwith open('static/basket_no2_2022.json') as f:\n  basket = json.load(f)\n\n# get the list of the urls in the NO2 files in netCDF format\ndf = pd.DataFrame(basket).T\nfile_urls = df['opendapUrl'].tolist()\nprint(file_urls)\n

['https://thredds.nilu.no/thredds/dodsC/ebas_doi/JE/CK/XJ/JECK-XJ9C.nc', 'https://thredds.nilu.no/thredds/dodsC/ebas_doi/XC/VP/A5/XCVP-A5BH.nc', 'https://thredds.nilu.no/thredds/dodsC/ebas_doi/3M/9Z/64/3M9Z-648F.nc', 'https://thredds.nilu.no/thredds/dodsC/ebas_doi/9K/6S/D6/9K6S-D6A9.nc', 'https://thredds.nilu.no/thredds/dodsC/ebas_doi/DQ/83/EQ/DQ83-EQBE.nc']\n

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

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

Define the year and variable of interest\u00b6

In\u00a0[4]:

# define the variable of interest and the year of interest\nvariable_interested_standard_name = 'mass_concentration_of_nitrogen_dioxide_expressed_as_nitrogen_in_air'\nyear_interested = 2022\n

Calculate daily, monthly and annual mean concentrations\u00b6

Inspect the data coverage by printing out how many days of data are available

In\u00a0[5]:

# define dataframes with the data of interest\ndaily_means_all_sites = []\nmonthly_means_all_sites = []\nannual_means_all_sites = []\nsite_name_list = []\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's time range overlaps with the wanted period\n    dataset_start_time = data.time.min().values\n    dataset_end_time = data.time.max().values\n\n\n    # check time period of the data\n    if (data.time.dt.year == year_interested).any():\n\n        # check if the dataset contains the variable of interest\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 the year of interest\n                variable_values = data[var].sel(time = data.time.dt.year == year_interested)\n\n                # Calculate daily, monthly and annual means\n                daily_means = variable_values.resample(time='1D').mean('time')\n                print('      ..has data for '+str(len(daily_means))+' days..')\n                monthly_means = variable_values.resample(time='ME').mean('time')\n                annual_means = variable_values.resample(time='YE').mean('time')\n\n                #add the means to the file that collects all the sites\n                daily_means_all_sites.append(daily_means)\n                monthly_means_all_sites.append(monthly_means)\n                annual_means_all_sites.append(annual_means)\n                site_name_list.append(data.attrs['ebas_station_name'])\n                break              \n\n        else:\n            print(f\"File {i} does not contain data in the desired time range.\")\n    data.close()\n

0\n File 0, FI0050R, variable: v_nitrogen_dioxide_ug_N_per_m3_amean (Hyyti\u00e4l\u00e4)                               Period: 2018-01-01T00:30:00.000000000 - 2022-12-31T23:30:00.000000000\n      ..has data for 365 days..\n1\n File 1, IT0009R, variable: v_nitrogen_dioxide_ug_N_per_m3_amean (Monte Cimone)                           Period: 2021-01-01T00:30:00.000000000 - 2022-12-31T23:30:00.000000000\n      ..has data for 365 days..\n2\n File 2, CZ0003R, variable: v_nitrogen_dioxide_ug_N_per_m3      (Kosetice (NAOK))                        Period: 2022-01-01T00:30:00.000000000 - 2022-12-31T23:30:00.000000000\n      ..has data for 365 days..\n3\n File 3, FR0020R, variable: v_nitrogen_dioxide_ug_N_per_m3      (SIRTA Atmospheric Research Observatory) Period: 2022-01-01T00:30:00.000000000 - 2022-12-31T23:30:00.000000000\n      ..has data for 365 days..\n4\n File 4, CH0001G, variable: v_nitrogen_dioxide_ug_N_per_m3      (Jungfraujoch)                           Period: 2021-01-01T00:30:00.000000000 - 2022-12-31T23:30:00.000000000\n      ..has data for 365 days..\n

Plot daily and monthly with one panel for each site\u00b6

In\u00a0[6]:

#First, prepare some variables to use as x values in the plotting\nndays=len(daily_means.time)\nnmonths=len(monthly_means.time)\ndays=np.linspace(1,ndays,ndays)\nmonthsdaily=np.linspace(1,nmonths,ndays)    #x values for the daily mean dots in the plot\nmonths=np.linspace(1,nmonths,nmonths)       #x values for the monthly mean lines in the plot\nnstations=len(site_name_list)\n\n\n# Prepare the figure. panelheight can be varied.\npanelheight=3\nfig, axes = plt.subplots(nrows=nstations,ncols=1, figsize=(10,int(panelheight*len(site_name_list))))\n\nfor col in range(0,nstations):\n    ax=axes[col]\n    ax.scatter(monthsdaily,daily_means_all_sites[col],label='Daily means, '+str(year_interested),alpha=0.3) #scatterplot of daily mean values\n    ax.plot(months,monthly_means_all_sites[col],label='Monthly means, '+str(year_interested),linewidth=2) #line plot of monthly mean values\n    ax.set_ylabel('NO$_2$ conc. (\u00b5gN/m$^3$)')\n    ax.set_title(site_name_list[col])\n    ax.legend(loc='upper left')\n    ax.set_xticks(np.arange(1, 13), ['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec'])\n\n\nplt.tight_layout()\nplt.show()\n

\n $\"No$ \n

Plot the monthly mean data from all the sites in the same panel\u00b6

In\u00a0[7]:

plt.figure(figsize=(10, 6))\nfor i in range(len(monthly_means_all_sites)):\n    plt.plot(monthly_means_all_sites[i], label=site_name_list[i], marker='*')\n\nplt.legend()\nplt.ylabel('\u00b5gN/m$^3$')\nplt.title('Monthly mean concentrations of NO$_2$  for year %d' % year_interested)\nplt.xticks(np.arange(0, 12), ['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec'])\nplt.show()\n

\n $\"No$ \n

Make a table with the aggregated data\u00b6

In\u00a0[8]:

#Define a file to included in the table\ndata = []\n\nMonths = ['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec']\n\nfor site_monthly, site_annual in zip(monthly_means_all_sites, annual_means_all_sites):\n    site_name = site_monthly.attrs['ebas_station_name']\n    lat = round(float(site_monthly.attrs['ebas_station_latitude']), 2)\n    lon = round(float(site_monthly.attrs['ebas_station_longitude']),2)\n    monthly_means = [round(mean, 2) for mean in site_monthly.values.tolist()] # Extract and round monthly means\n    annual_mean = [round(mean, 2) for mean in site_annual.values.tolist()] # Extract and round annual mean\n     # Append data to the list\n    data.append([site_name, lat, lon] + monthly_means + annual_mean) \n\n# Define column names\ncolumns = ['Site Name', 'Latitude', 'Longitude'] + Months + ['Annual Mean']\n\n# Create DataFrame\ndf_no2_2022 = pd.DataFrame(data, columns=columns)\ndf_no2_2022\n

Out[8]:

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

	Site Name	Latitude	Longitude	Jan	Feb	Mar	Apr	May	Jun	Jul	Aug	Sep	Oct	Nov	Dec	Annual Mean
0	Hyyti\u00e4l\u00e4	61.85	24.28	0.60	0.79	0.48	0.21	0.19	0.28	0.24	0.34	0.31	0.45	0.75	0.95	0.47
1	Monte Cimone	44.19	10.70	0.19	0.28	0.23	0.20	0.12	0.12	0.14	0.20	0.15	0.08	0.19	0.16	0.17
2	Kosetice (NAOK)	49.57	15.08	4.22	2.83	4.34	3.52	3.01	2.85	2.85	3.34	3.17	3.38	3.93	6.23	3.62
3	SIRTA Atmospheric Research Observatory	48.71	2.15	3.92	1.43	4.17	2.11	1.49	1.82	1.81	2.01	2.40	2.30	2.12	4.49	2.53
4	Jungfraujoch	46.55	7.99	0.04	0.06	0.06	0.11	0.07	0.06	0.06	0.07	0.06	0.05	0.05	NaN	0.06

Make a map plot with the annual means colored by its average value\u00b6

In\u00a0[9]:

#Extract data from the df_no2_2022 DataFrame created above\nsite_name_list = df_no2_2022['Site Name'].tolist()\nlongitudes = df_no2_2022['Longitude'].tolist()\nlatitudes = df_no2_2022['Latitude'].tolist()\nmean_values = df_no2_2022['Annual Mean'].tolist()\n\n#Choose max and min - this will depend on your data\nminval=0\nmaxval=5 \n\n# Create the map plot\ntiler = cimgt.OSM()\nfig = plt.figure(figsize=(9, 9))\nax = plt.axes(projection=tiler.crs)\n\n# Add map features (coastlines, borders, etc.)\nzoom = 4\nax.add_image(tiler, zoom)\n\n# Define discrete color levels and the colormap\ncmap = cm.get_cmap('viridis').reversed()  # Choose your colormap, e.g., 'viridis', 'plasma', etc.\nncolors = 15  # 15 color levels\nlevels = np.linspace(minval, maxval, ncolors)  # Discrete values for the colorbar (e.g., from 0 to 1)\nnorm = BoundaryNorm(boundaries=levels, ncolors=cmap.N)\n\n# Plot the stations as colored circles\nsc = ax.scatter(longitudes, latitudes, c=mean_values, cmap=cmap, norm=norm, s=70, edgecolor='black', transform=ccrs.PlateCarree())\n\n# Numerate the stations on the map\nfor i, (lon, lat) in enumerate(zip(longitudes, latitudes), start=1):\n    ax.text(lon + 0.8, lat + 0.8, str(i), fontsize=16, ha='left', va='center', \n            transform=ccrs.PlateCarree())  # , fontweight='bold')\n\n# Add a colorbar with discrete intervals\ncbar = plt.colorbar(sc, ax=ax, orientation='vertical', pad=0.05)\ncbar.set_label('\u00b5gN/m$^3$')\ncbar.set_ticks(levels)  # Set colorbar ticks to discrete levels\ncbar.ax.set_yticklabels([f'{level:.1f}' for level in levels])  # Customize tick labels if needed\n\n# Set titles, labels, and other map properties\nax.set_title('NO$_2$ annual mean concentrations in ' + str(year_interested))\n\nax.set_extent([min(longitudes) - 10, max(longitudes) + 10, min(latitudes) - 10, max(latitudes) + 10], crs=ccrs.PlateCarree())  # Set the extent of the map\n\n# Add the legend as text, positioned just below the map\nlegend_text = \", \".join([f\"{i}. {name}\" for i, name in enumerate(site_name_list, start=1)])\nplt.figtext(0.1, 0.01, legend_text, ha=\"left\", fontsize=9, wrap=True)\n\n# Adjust the layout to avoid cutting off the legend\nplt.tight_layout(rect=[0.01, 0.01, 1, 1])\nplt.show()\n

\n $\"No$ \n

To further analyse NO2 data by looking at trends and climatology, have a look at notebook example \"Trends and clomatology NO2.ipynb\".