FZ filter not running

[MarisZul]

The FZ filter seems to work on small datasets such as Amsterdam and Götheborg dataset. For two other, larger datasets the following behaviour was observed:

Dataset 1 in hourly resolution gives the following error:

Strangely, the same Data in 5 min resolution works.

Dataset 2 in hourly resolution gives the same error.

It seems that after adding the private function of nbrs_not_nan the issue arose. When resolving the issue pay attention to:
1.) the private function uses nr_nbrs_not_nan = (ds_pws.rainfall.sel(id=neighbor_ids).**notnull**().sum(dim="id")) while the old version uses **.isnull()**
2.) the new private function uses the length of the time dim to create an array, which might in some cases cause the "index out of bounds problem"

[MarisZul]

When the same Dataset is transposed such that the variable is rainfall(id, time) the filter works! Just a reminder that when using the FZ filter the dimensions have to follow the right order. This should be included in the notebook and maybe the function could be adapted to be more flexible i.e allowing both possibilities.

[cchwala]

Thanks for figuring this out! IMO it would be best to always transpose the Dataset or DataArray inside the function since we rely on fixed naming of the dimension anyway we can also enforce a certain order.

[MarisZul]

Good idea! I will add another private function that checks if the dimension order is correct and otherwise transposes to the desired order. The private function will be called in every function. Then I will write a test for the private function with an example dataset with "wrong" dimension order to achieve full patch coverage. I plan to to this by the end of the week with another PR.
Let me know, if you have any other ideas/suggestions.

[cchwala]

...that checks if the dimension order is correct and otherwise transposes to the desired order.

With transpose in xarray you can provide an order of named dimension. Hence, you do not need to first check the order of dims and then decided if to transpose or not. You can just always call this function with the required order of dims. I assume this function is computationally very cheap since it will not transpose any data, just move the reference to the dims.

[MarisZul]

You´re right, the transpose will do nothing if the order is already correct! Then adding one line ds_pws = ds_pws.transpose("id","time") to each function is enough and it will run either way so no tests neccessary.

[cchwala]

It might still be a good idea to have a function that checks the dataset. Actually, we should have done that already some time ago, maybe in poligrain because, if I recall correctly, some people already had similar issue with wrong dim order.

Hence, what you could do is to add a function that does something like a "structure check" (but find a better name) of the dataset. You could then call this function everywhere where it is needed, probably at the top of all the high-level functions that take a dataset as input where the order of dims matter, or maybe for all functions that take a dataset as input because you never know if the dim order will matter later on.

[cchwala]

With "structure check" I meant something like I describe in OpenSenseAction/poligrain#38. But if it is an all-encompassing data format check, this is definitely a functionality that has to be implemented in poligrain and not in pypwsqc.

I am not sure how to do this best. I hesitate a bit with the idea to call a function to check the structure of a Dataset in every function that does some processing. I will look at the code now and see if I have an idea or preference.

[cchwala]

The code that causes the problems is this one

pypwsqc/src/pypwsqc/flagging.py

Lines 71 to 85 in ca2025c

	for i in np.arange(len(pws_data.id.data)):
	for j in np.arange(len(pws_data.time.data)):
	if j < nint:
	fz_array[i, j] = -1
	elif sensor_array[i, j] > 0:
	fz_array[i, j] = 0
	elif fz_array[i, j - 1] == 1:
	fz_array[i, j] = 1
	elif (np.sum(sensor_array[i, j - nint : j + 1]) > 0) or (
	np.sum(ref_array[i, j - nint : j + 1]) < nint + 1
	):
	fz_array[i, j] = 0
	else:
	fz_array[i, j] = 1

This code is special because it is the only part in flagging.py that use numpy-arrays and direct 2D-indexing of these where the order of dims matters. Based on a quick check of the other modules in pypwsqc I would say this is currently the only part of the code which is sensitive to the dim order.

Hence, I would just fix it here by adding a transpose using xarray at the top of the function.

Looking again at this nested for-loop, I would say, this part of the code could be optimized to a vectorized form without for loop.

If and how a structure check of a Dataset is done, we should discuss and decide in OpenSenseAction/poligrain#38.