pyaerocom icon indicating copy to clipboard operation
pyaerocom copied to clipboard

MDA8

Open WillemVanCaspel opened this issue 2 years ago • 8 comments

Hi everyone,

The daily maximum of the 8-hourly running mean ozone concentration (MDA8) is a metric often used by the WHO, and it would be really great if the observed values could be compared to EMEP (the model already outputs it). Is there any chance it would be possible to evaluate MDA8 using pyaerocom?

Willem

WillemVanCaspel avatar May 22 '23 08:05 WillemVanCaspel

In principle this should be easy to add since we calculate the variable vmro3max e.g. at EBAS (grepping through the code for vmro3max reveals this):

data/variables.ini:[vmro3max]
data/variables.ini:var_name = vmro3max
data/emep_variables.ini:vmro3max = "SURF_MAXO3"
data/ebas_config.ini:[vmro3max]
aux_var_helpers.py:def calc_vmro3max(data):
aux_var_helpers.py:    new_var_name = "vmro3max"
aux_var_helpers.py:    o3max = data[var_name]
aux_var_helpers.py:    return o3max
plugins/mep/aux_vars.py:def vmro3max_from_ds(ds: xr.Dataset) -> xr.DataArray:
plugins/mep/aux_vars.py:    return conc_to_vmr(ds["conco3"], vmr="vmro3max")
plugins/mep/reader.py:from .aux_vars import vmrno2_from_ds, vmro3_from_ds, vmro3max_from_ds
plugins/mep/reader.py:        "vmro3max": vmro3max_from_ds,
plugins/mscw_ctm/emep_variables.toml:vmro3max = "SURF_MAXO3"
aeroval/glob_defaults.py:    "vmro3max": {"scale": [0, 7.5, 15, 22.5, 30, 37.5, 45, 52.5, 60], "colmap": "coolwarm"},
aeroval/glob_defaults.py:    vmro3max=["O3Max", "3D", "Volume mixing ratios"],
io/read_ebas.py:    calc_vmro3max,
io/read_ebas.py:        "vmro3max": ["vmro3"],
io/read_ebas.py:        "vmro3max": calc_vmro3max,
io/read_eea_aqerep_base.py:    VAR_NAMES_FILE["vmro3max"] = "concentration"
io/read_eea_aqerep_base.py:        vmro3max="**/??_7_*_timeseries.csv*",
io/read_eea_aqerep_base.py:    CONV_FACTOR["vmro3max"] = np.float_(
io/read_eea_aqerep_base.py:    CONV_UNIT["vmro3max"] = "ppb"
io/read_eea_aqerep_base.py:        "vmro3max": ["conco3"],
io/read_eea_aqerep_base.py:        "vmro3max": NotImplementedError(),

so adding the calculation of the 8hourly running mean should be easy.

I suggest vmrmda8 as variable name

jgriesfeller avatar May 23 '23 11:05 jgriesfeller

Check the WHO definition for how this is computed. What is the minimum temporal coverage required to compute this statistics?

lewisblake avatar Jun 19 '23 09:06 lewisblake

Which project should allocate resources for this / how badly is it needed?

lewisblake avatar Oct 23 '23 09:10 lewisblake

It's not needed very badly. Since posting this issue, I have made my own post-processing script that calculates peak season MDA8 from the co-location data files produced by Pyaerocom.

WillemVanCaspel avatar Oct 23 '23 10:10 WillemVanCaspel

Michael G says:

yes, MDA8 is needed for the FAIRMODE plots

I can also confirm that the 8-hour average is only to be calculated if at least 6 out of 8 measurements are available

what I haven't found out yet is which 8-hour means are to be calculated during a day. There are different methods found on the internet. E.g. 0-7, 1-8, ... 16-23, i.e. 17 averages per day. Some authors, as we know, consider a day to last from 6 UTC to 6 UTC. I've also read somewhere that there are 24 8-hour averages during one day (which sounds weird..) oh dear... I'll ask Meteo France what they've done during the last 10 years, and we'll do the same for continuity. I'll let you know.

lewisblake avatar Feb 07 '24 12:02 lewisblake

From the EU directive document annex XI:

Eight hours values:         75 % of values (i.e. 6 hours)

Maximum daily 8-hour mean:  75 % of the hourly running eight hour
                            averages (i.e. 18 eight hour averages per
                            day)
                            
The maximum daily eight hour mean concentration will be selected by examining 
eight hour running averages, calculated from hourly data and updated each hour. 
Each eight hour average so calculated will be assigned to the day on which it 
ends i.e. the first calculation period for any one day will be the period from 
17:00 on the previous day to 01:00 on that day; the last calculation period for 
any one day will be the period from 16:00 to 24:00 on that day.

WillemVanCaspel avatar Feb 07 '24 12:02 WillemVanCaspel

yeah, sorry I was slow, had a visitor over lunch. So Willem is write, and this is also what MeteoF is doing.

michaelgau avatar Feb 07 '24 13:02 michaelgau

In CAMS2_83 we have to evaluate 1-day analyses and 4-day forecasts. The 8-hour averages starting on the previous day cannot be processed in the case of the first days. Thus I think for the one day of the analysis and the 1st day of the forecast we can calculate only 17 averages. ...unless we mix data from analyses/forecasts that are issued on different days, but that would make it a bit messy.

michaelgau avatar Feb 07 '24 13:02 michaelgau