Issue on some weeks and points in observational data: hindcast-like-observations_2000-2019_biweekly_terciled.zarr

Hi! I'm working on the challenge and I have one question about the observation data. After plotting the data for the day 02-27 over all years I can see that there are some points where there a lot of observations corresponding to one category, for example the point with latitude=25.5 and longitude=43.5 only have two observations on the category near normal. This fact happens on many other points and also in others days of the year. Aren't observations supposed to be better distributed? I would expect to see a more uniform distribution around the three categories.

I attach images of the code that I use and the graphics that illustrate this data.

Edit: If it is useful for someone the same strange thing occur if you check the observation data for 10 Dec. There's a point with latitude=60 longitude=13.5 that have no near normal observations along all 20 years.

first look: https://gist.github.com/aaronspring/60ff16e58c1debd2e88f38cc7f265cea

My explanation: category edges are calculated from observation for each week, see groupby('forecast_time.weekofyear') in https://renkulab.io/gitlab/aaron.spring/s2s-ai-challenge-template/-/blob/master/notebooks/renku_datasets_biweekly.ipynb.

why each week? I didnt choose to calculate by dayofyear, because in leap years dayofyear is shifted by +1 starting in march.

changed the description

maybe we should check the original raw observations before categorization and the edges

https://gist.github.com/aaronspring/60ff16e58c1debd2e88f38cc7f265cea

Some more context https://github.com/pydata/xarray/issues/1844

So the main issue here, as I see it now, is the leap day in the years 2004, 2008, 2012, 2016 and 2020, which shifts forecast_time.dayofyear by +1 and also shifts the forecast_time.dt.weekofyear by +1. I had hoped to circumvent that issue by taking weekofyear, but didnt pay enough attention.

Solution ideas:

• One practical way is to just drop Feb 29th in S2S. However, then we are missing on day or have a gap. Not an option here IMO.
• The observations clearly contain Feb 29th. Another way of dealing with this is triangular smoothing in the climatology, i.e. taking more than just dayofyear or weekofyear into account.
• Just ignore it. Most grid cells dont show too large differences in the freq distributions between categories.

Hi @aaron.spring this is Llorenç from BSC, working with Sergi. The problem seems to be that weekofyear uses the ISO standard, which says that the first week of the year is the first week with a Thursday on that year. This means that 2nd January is sometimes 1st week, or other times 52nd week or even 53rd week, depending on what day is January 1st:

Droping 29th Feb will not solve the issue, I believe. If you are interested, we have some code to split forecast_time into (year, week) components, which can be used later to group by.

We believe this is a very important issue to solve. Having categories that are not equiprobable makes the forecasting problem much difficult for "classical methods" that count number of members.

Thanks for the example for Jan 2nd also @llledo More context: https://en.wikipedia.org/wiki/ISO_week_date.

We believe this is a very important issue to solve. Having categories that are not equiprobable makes the forecasting problem much difficult for "classical methods" that count number of members.

I agree.

One way to go would be using dayofyear, but then we would have the 365/366 day issue. Another way would be splitting weeks manually by hand and then use this as a groupby, as you are suggesting. I'd be keen on seeing how you would split and groupby.

My first idea would be to have 365 days in the year and have the leap week end of feb having 8 days in leap years, similar to https://github.com/pydata/xarray/issues/1844.

Hi @aaron.spring here is the code we have to divide the start date into year and week components and then compute the tercile edges.

import xarray as xr

import numpy as np

import pandas as pd

obs_2000_2019 = xr.open_zarr("/path/to/hindcast-like-observations_2000-2019_biweekly_deterministic.zarr", consolidated=True)

y = list(range(2000,2019+1))

w = list(range(1,54))

tmp = obs_2000_2019.forecast_time.values.reshape(20,53)

obs_2000_2019 = obs_2000_2019.assign_coords(year=y, week=w).stack(forecast_time=("year", "week")).unstack("forecast_time")

obs_2000_2019.coords["forecast_time"] = (["year", "week"], tmp)

edges = obs_2000_2019.chunk({'year':-1,'week':-1,'longitude':'auto'}).groupby('week').quantile(q=[1./3.,2./3.], dim=['year']).rename({'quantile':'category_edge'}).astype('float32')

obs_tercile_2000_2019 = make_probabilistic(obs_2000_2019, edges)

Thank you very much for this code snippet @sbechsal7. I will mention the two of you in the next release, where these categorizations will be fixed.

So basically you are breaking 1d-forecast_time into 2d-forecast_time(year, week). I will have to think how my notebooks will still work downstream and I'm confident this can be easily integrated.

Can you confirm that the category distribution is now more equalized? I guess so - counting 6 or 7 occurrences in 20 years per category?

Yes, we made some plots and it seems that the category distribution is equalized.

mentioned in merge request s2s-ai-challenge-template!18 (merged)

assigned to @aaron.spring

closed with merge request s2s-ai-challenge-template!18 (merged)

mentioned in commit s2s-ai-challenge-template@40fc39be