![\"No](\"https://raw.githubusercontent.com/ACTRIS-Data-Centre/actris-jupyter-hub/refs/heads/main/actris-dc-logo.png\"){kind=logo}
\nThe ATMO-ACCESS Virtual Access Portal gives access to Footprint analysis tools for ground based aerosol measurements. Through the portal you can request model runs to produce data products (e.g. footprint, source contributions) for your decided locations, or search and use the already produced products.\nAccess the ATMO-ACCESS Virtual Services: https://www.atmo-access.eu/virtual-access/
\nThis example plots the footprint emission sensitivity for analysing black carbon (BC) measurements at a given facility. The products are produced using FLEXPART model and include footprint analysis for the air mass arriving at the station on 3 hour time resolution 20 days back in time. The footprint emissions sensitivity showing the probability of any release occurring in any grid-cell to reach the receptor (station) for 30 days particle tracking. There is also modelled source analysis available for some of the products. More documentation of the products are available at https://flexpart-request.nilu.no/data-access.
\nIf you are interested in combining timeseries of modelled source analysis of black carbon with measured black carbon, have a look at example \"FLEXPART and observed timeseries.ipynb\"
\nImport libraries\u00b6
\nimport sys\nimport numpy as np\nfrom netCDF4 import Dataset\nimport requests\nimport io\nfrom datetime import datetime, timedelta\nimport matplotlib.pyplot as plt\nfrom mpl_toolkits.basemap import Basemap, cm, maskoceans\nfrom matplotlib.colors import LinearSegmentedColormap\n/home/lemu/miniconda3/envs/actris-vre/lib/python3.14/site-packages/pyproj/network.py:59: UserWarning: pyproj unable to set PROJ database path.\n _set_context_ca_bundle_path(ca_bundle_path)\n\n
Settings for FLEXPART footprint emission sensitvity\u00b6
\nChoose the settings you want for your FLEXPART footprint emission sensitivity plot. You can choose facility, year, month, altidute, day and hourly interval. Lastly you can choose the projection you want for the footprint.
\nFacilities to choose from:
\nNO0042G - Zeppelin mountain (Ny-\u00c5lesund) *(https://data.actris.eu/facility/w2kl)*
\nFI0050R - Hyyti\u00e4l\u00e4 *(https://data.actris.eu/facility/ko6m)*
\nFR0030R - Puy de D\u00f4me *(https://data.actris.eu/facility/do7b)*
\nCZ0003R - Kosetice (NAOK) *(https://data.actris.eu/facility/p797)*
\nGR0002R - Finokalia *(https://data.actris.eu/facility/wdd9)*
\nCH0002R - Payerne *(https://data.actris.eu/facility/7tbf)*
\nCH0001G - Jungfraujoch *(https://data.actris.eu/facility/mmee)*
\nGR0100B - Athens Demokritos *(https://data.actris.eu/facility/ekme)*
\nES0018G - Izana *(https://data.actris.eu/facility/flqs)*
\nIT0009R - Monte Cimone *(https://data.actris.eu/facility/7uqv)*
\nNO0002R - Birkenes II *(https://data.actris.eu/facility/9cxe)*
\nES0100R - El Arenosillo *(https://data.actris.eu/facility/3pb5)*
\nGR0101R - Helmos Mountain *(https://data.actris.eu/facility/n1rx)*
\nDE0043G - Hohenpeissenberg *(https://data.actris.eu/facility/rhhz)*
\nDE0044R - Melpitz *(https://data.actris.eu/facility/fhhq)*
\nES0019U - Barcelona (Palau Reial) *(https://data.actris.eu/facility/hcwo)*
\nES0020U - Granada *(https://data.actris.eu/facility/5qss)*
\nFR0020R - SIRTA - Orme des Merisiers *(https://data.actris.eu/facility/aswe)*
\nIT0004R - Ispra *(https://data.actris.eu/facility/j133)*
\nIT0013R - Potenza *(https://data.actris.eu/facility/epdw)*
\nRO0007R - Bucharest I *(https://data.actris.eu/facility/ab2x)*
\nFR0026R - La R\u00e9union - Ma\u00efdo atmospheric observatory *(https://data.actris.eu/facility/ka8v)*
\n# SETTINGS, FEEL FREE TO CHANGE AND ADAPT TO YOUR NEEDS\nSTATION = \"DE0043G\" # GIVE STATION CODE\n# NOTE stations DE0043G\u2006 \u2006 DE0044R\u2006 \u2006 ES0019U\u2006 \u2006 ES0020U\u2006 \u2006 FR0020R\u2006 \u2006 IT0004R\u2006 \u2006 IT0013R\u2006 \u2006 RO0007R FR0026R are lidar and for these ALT>0 is allowed [0, 100, 500, 4000]\nYEAR = 2023# GIVE YEAR\nMON = 7 # GIVE MONTH\nALT = 100 # GIVE ALTITUDE\nDD = 20 # GIVE DAY\nTT = 12 # GIVE THE TIME IN 3 HOURLY INTERVALS LIKE 0, 3, 6, ..., 21\nPROJ = \"REGIONAL\" # GIVE THE PROJECTION YOU WANT TO SEE: CHOOSE BETWEEN GLOBAL, REGIONAL, POLAR \n\n##### DO NOT CHANGE BELOW THIS UNLESS YOU KNOW WHAT YOU DO #####\n# Available stations and their coordinates\nSTNLIST = [\"NO0042G\", \"FI0050R\", \"FR0030R\", \"CZ0003R\", \"GR0002R\", \"CH0002R\", \"CH0001G\", \"GR0100B\", \"ES0018G\", \"IT0009R\", \"NO0002R\", \"ES0100R\", \"GR0101R\",\"DE0043G\", \"DE0044R\", \"ES0019U\", \"ES0020U\", \"FR0020R\", \"IT0004R\", \"IT0013R\", \"RO0007R\", \"FR0026R\"]\nLONLIST = [11.88668, 24.283333, 2.964886, 15.080278, 25.669399, 6.944722, 7.985, 23.816, -16.4994, 10.7, 8.252, -6.733, 22.196262, 11.009619, 12.9277, 2.116111, -3.605, 2.158889, 8.633333, 15.7237, 26.029, 55.383]\nLATLIST = [78.90715, 61.85, 45.772223, 49.573394, 35.337799, 46.813056, 46.5475, 37.994999, 28.309, 44.183333, 58.38853, 37.1, 37.984265, 47.801498, 51.5255, 41.39, 37.164, 48.708611, 45.8, 40.601, 44.348, -21.0794]\nALTLIST = [0, 100, 500, 4000, 6000]\n# Find coordinates of the station given in SETTINGS\nxx = []; yy = []\nfor i in range(len(STNLIST)):\n if STNLIST[i] == STATION:\n xx.append(LONLIST[i])\n yy.append(LATLIST[i])\n\nprint(\"The coordinates of\",STATION, \"are\",xx, \"and\",yy)\nThe coordinates of DE0043G are [11.009619] and [47.801498]\n\n
# Find altitude of the station given in SETTINGS\nzz1 = []; zz2 = []\nfor i in range(len(ALTLIST)):\n if ALTLIST[i] == ALT:\n zz1.append(ALTLIST[i])\n zz2.append(ALTLIST[i+1])\n\nprint(\"The altitude of\",STATION, \"is\",zz1, \"and\",zz2)\nThe altitude of DE0043G is [100] and [500]\n\n
# Set datetime objects dat and datE\ndat = datetime(YEAR,MON,DD,TT)\ndatE = dat + timedelta(hours=3)\n\n# Define the start and end of the month\nstart_date = datetime(YEAR, MON, 1, 0, 0) # Start at midnight on the first day\n# Find the first day of the next month to mark the end\nif MON == 12:\n end_date = datetime(YEAR + 1, 1, 1, 0, 0) # If December, next month is January of the next year\nelse:\n end_date = datetime(YEAR, MON + 1, 1, 0, 0) # Otherwise, just increment the month\n\ntime_steps = []\ncurrent_time = start_date\nwhile current_time < end_date:\n time_steps.append(current_time)\n current_time += timedelta(hours=3)\n\n\n#print(ddlist)\ntarget_datetime = datetime(YEAR,MON,DD,TT)\n# Find the index of the specific datetime within the list\ntry:\n idx = time_steps.index(target_datetime)+1\n print(f\"The index of {target_datetime} is {idx}.\")\nexcept ValueError:\n print(f\"{target_datetime} is not in the list of time steps.\")\n\nprint(\"The index of \", dat, \"within a month is\", idx)\nThe index of 2023-07-20 12:00:00 is 157.\nThe index of 2023-07-20 12:00:00 within a month is 157\n\n
# URL of the NetCDF file (example URL)\nurl = 'https://atmo-access.nilu.no/WEBPLOTS/' + str(STATION) + '/' + str(YEAR) + '/' + \\\n'OUT_' + str(YEAR) + str(MON).zfill(2) + '_' + str(ALT) + '/' + 'ECMWF/regional_footprint_a/OUT_'+ \\\nstr(YEAR) + str(MON).zfill(2) + '_' + str(ALT) +'.regional_footprint_a_'+str(idx)+'.nc'\nprint(\"File is taken from\", url)\nFile is taken from https://atmo-access.nilu.no/WEBPLOTS/DE0043G/2023/OUT_202307_100/ECMWF/regional_footprint_a/OUT_202307_100.regional_footprint_a_157.nc\n\n
# Function to read the NetCDF file from the internet\ndef read_netcdf_from_url(url):\n try:\n # Send a GET request to download the file\n response = requests.get(url)\n response.raise_for_status() # Raise an exception for HTTP errors\n \n # Open the file from the response content using an in-memory bytes buffer\n with io.BytesIO(response.content) as data:\n # Open the dataset from the byte stream\n ds = Dataset('inmemory.nc', mode='r', memory=data.read())\n # Print basic information about the dataset\n print(ds) \n # Example: List variables in the dataset\n print(\"Variables in the dataset:\")\n print(list(ds.variables.keys()))\n lons = ds.variables['longitude'][:]\n lats = ds.variables['latitude'][:]\n foot = ds.variables['footprint'][:]\n \n # Close the dataset when done\n ds.close()\n print(\"NetCDF file processed successfully.\")\n except requests.exceptions.RequestException as e:\n print(f\"Error downloading the NetCDF file: {e}\")\n except Exception as e:\n print(f\"Error reading the NetCDF file: {e}\")\n return lons, lats, foot\n\n# Call the function to read the NetCDF file from the URL\nlons, lats, foot = read_netcdf_from_url(url)\n<class 'netCDF4.Dataset'>\nroot group (NETCDF4 data model, file format HDF5):\n author: Nikolaos Evangeliou\n contact: ne@nilu.no\n created_date_UTC: 19-03-2024\n created_hour_UTC: 02:49:02\n institution: Norsk Institute for Luftforskning\n description: Creates 2D netCDF from a pre-existing array\n history: Created by NE\n source: netCDF4 python module\n coordinates: latitude longitude\n dimensions(sizes): latitude(360), longitude(720)\n variables(dimensions): float32 latitude(latitude), float32 longitude(longitude), float32 footprint(latitude, longitude)\n groups: \nVariables in the dataset:\n['latitude', 'longitude', 'footprint']\nNetCDF file processed successfully.\n\n
######## MAKE PLOT ########\ndef plot(dat,datE,xx,yy,lons,lats,zpoint1,zpoint2,var,area='GLOBAL'):\n fig = plt.figure(figsize=(8,7.5)) # Ensure fig is defined\n ax = fig.add_subplot(111)\n ax.set_title('Footprint emission sensitivity for BC'+'\\nfrom '+ dat.strftime(\"%d-%b-%Y %H:%M:%S\")+' to '+datE.strftime(\"%d-%b-%Y %H:%M:%S\")+'\\nLowest release height: '+str(int(zpoint1))+' m Highest release height: '+str(int(zpoint2))+' m',fontsize=14)\n \n if area == 'GLOBAL':\n m = Basemap(projection='robin',lon_0=0,resolution='l')\n\n elif area == 'REGIONAL':\n m = Basemap(width=8000000,height=5000000,resolution='l',projection='stere',lat_ts=60.,lat_0=yy,lon_0=xx)\n\n elif area == 'POLAR':\n h = 10000.\n m = Basemap(projection='nsper',lon_0=xx,lat_0=90.,satellite_height=h*1000.,resolution='l')\n\n else:\n print(\"Error\")\n sys.exit()\n\n lon, lat = np.meshgrid(lons, lats)\n xi, yi = m(lon, lat)\n \n levs = np.arange(-2,6.01,0.01)\n clevs = 2**levs\n boundlevs = [str(('%.2f' % v).rstrip('0').rstrip('.')) if v <=1 else str(int(round(v))) for v in clevs]\n \n m.drawcoastlines()\n if area == 'GLOBAL' or area == 'POLAR':\n m.drawparallels(np.arange(-90.,91.,30.), labels=[1,0,0,0], fontsize=10)\n m.drawmeridians(np.arange(-180.,181.,60.), labels=[0,0,0,1], fontsize=10)\n elif area == 'REGIONAL':\n m.drawparallels(np.arange(-90.,91.,10.), labels=[1,0,0,0], fontsize=10)\n m.drawmeridians(np.arange(-180.,181.,20.), labels=[0,0,0,1], fontsize=10)\n m.drawmapboundary()\n m.drawcountries()\n \n colors = [\"purple\", \"red\", \"orange\", \"gold\", \"lawngreen\", \"cyan\", \"blue\", \"mediumorchid\", \"white\"]\n cmap = LinearSegmentedColormap.from_list(\"cb_friendly\", colors[::-1], N=900)\n norm = cm.colors.BoundaryNorm(boundaries=clevs,ncolors=cmap.N)\n\n cs = ax.contourf(xi, yi, var[:, :], clevs, cmap=cmap, extend='both', norm=norm)\n cbar = m.colorbar(cs,ticks=clevs[::100],location='bottom', pad=\"12%\")\n cbar.set_label('seconds',fontsize=15)\n cbar.ax.set_xticklabels(boundlevs[::100])\n \n aa1,bb1= m(xx, yy)\n m.plot(aa1, bb1, marker='*', markersize=10, markeredgecolor='k',markeredgewidth=0.8,markerfacecolor='m',linestyle='')\n plt.show()\n\n# PLOT IT\u2006 \u2006 \u2006 \u2006 \nplot(dat,datE,xx[0],yy[0],lons,lats,zz1[0],zz2[0],foot,area=PROJ)\nIf you want to compare modelled source analysis of black carbon together with measured black carbon, have a look at example \"FLEXPART and observed timeseries.ipynb\"
\n