ESMValGroup / ESMValTool

ESMValTool: A community diagnostic and performance metrics tool for routine evaluation of Earth system models in CMIP
https://www.esmvaltool.org
Apache License 2.0
224 stars 128 forks source link

"non contiguous" coordinates in area-weighted resampling #3387

Open BSchilperoort opened 1 year ago

BSchilperoort commented 1 year ago

I am trying to update the hydrology/wflow diagnostic script for the Wflow.jl model, and encountered a bug/issue in regridding.

In this case I want to regrid data using the "area-weighted" method. The data comes from CMORized ERA5. However, when trying to use esmvaltool.preprocessing.regrid, it fails due to coordinates being non-contiguous.

It turns out that, because for my test case a shapefile-based cutout is used with longitude values [-5 -> 5]. This corresponds to 360 degree longitude values [0->5, 355->360].

This then fails in Iris here. As the coordinates are not contiguous in a 0-360 degree based system.

I.e. any (longitudinal) cutout that crosses the prime meridian will encounter this issue and cannot be regridded.

valeriupredoi commented 1 year ago

Hi @BSchilperoort cheers for opening this issue! Please see this issue https://github.com/ESMValGroup/ESMValCore/issues/2229 in which @rebeccaherman1 documents this type of behaviour as well. We will look into finding a workaround, but IMHO it's probably best handled upstream at iris - I'll get back to you once we've discussed the options :beer:

BSchilperoort commented 1 year ago

Thanks for the quick reply! I did not find that issue using my search terms but it's good to know I'm not the only one encountering this problem.

BSchilperoort commented 1 year ago

It turns out that by removing and guessing the bounds the error can be bypassed. I don't think it will necessarily result in a miscalculation but I'm not 100% sure.

cube.coord("longitude").bounds = None
cube.coord("longitude").guess_bounds()
BSchilperoort commented 1 year ago

A more robust fix I use now is to shift the coordinates using:

cube = cube.intersection(longitude=(-180., 180.))