Bbox fails when wrapping around dataset boundaries

[nilshempelmann]

Description

a bbox :

lon0 = -5 # Minimum longitude.
lon1 = 17.0 # Maximum longitude.
lat0 = 10.5 # Minimum latitude.
lat1 = 24.0 # Maximum latitude.

out = finch_i.subset_bbox(
    resource=ds[0], lon0=lon0, lon1=lon1, lat0=lat0, lat1=lat1)

failed with:

owslib.wps.WPSException : {'code': 'NoApplicableCode', 'locator': 'None', 'text': 'Process error: method=wps_xsubsetbbox.py._handler, line=170, msg=Input longitude bounds ([-5. 17.]) cross the 0 degree meridian but dataset longitudes are all positive.'}

[Zeitsperre]

This is more of a limitation of xclim itself. The way to work around this would be to either split the box into two boxes using positive-only coordinates or to make a polygon based on those coordinates that would be fed to the subset_shape-based process (which is unfortunately not very optimized).

[nilshempelmann]

can you include an argument 'evaluation' and appropriate handling in the beginning of the process which is doing the job for the user?

[Zeitsperre]

I'm not sure if I follow you exactly, what do you mean by 'evaluation'?

[nilshempelmann]

To evaluate the input setting. (The arguments or parameter)
And if there is a lon0 = -1 which will distract xclim that the process will handle it appropriately.

[Zeitsperre]

Clipping the data according to boundaries set with negative longitudes is an allowed case (e.g. lon0 = -80, lon1 = -73). The algorithm examines the boundaries of the data and will add the 360 degrees needed to make it consistent with the NetCDF. The issues arise when the bounding box crosses the min and max longitudes of the data (ie: 0 and 360, or -180,180 if the data is organized that way).

Like I mentioned online, this subsetting method is a bit hacky. Realigning all the longitudes of a netCDF is costly (0,360 to -180,180) and while this is handled in the subset_shape process, it's not at all efficient (relies on several workarounds using GeoPandas.

I can look into this today to see if I can do something similar, but the code is starting to get really spaghetti-like. I might need to add a new decorator or evaluation process to sort which subsetting algorithm to use. Not ideal.

[tomLandry]

Hello guys. We're starting to look at the sample subsetting code from xclim with @tlogan2000 and @davidcaron. I guess problems mentionned here for a bounding box will appear for polygonal subsetting. In the short term, I'm not sure we can solve the general case but let's keep in touch.

[Zeitsperre]

@tomLandry Think you could keep me in the loop on this? I wrote much of the handler logic so I can weigh in on the reasoning behind the decorators and the cases they handle.

[nilshempelmann]

@tomLandry right for polygon subset its necessary to be able to handle meridian and dayline as well
Good luck :-)

[Zeitsperre]

I hadn't even considered the dayline. That's a whole new can of worms. Definitely worth considering moving forward.

[tlogan2000]

For the meridian we should be able to apply the same logic as the subset.polygon ...

[tlogan2000]

To have something functionnal quickly we could simply probably call subset.poly on the bbox rectanlge when this error is raised

[huard]

How does ocgis handle this ?

[nilshempelmann]

We was running in similar issues with ocgis bird-house/flyingpigeon#126
no problems with meridian, but with the dayline.
Than @bekozi added a wrapper to be switched depending on which part of the world you would like to subset. (the bbox over the dayline or the other way round a bbox covering the rest of the world)
https://ocgis.readthedocs.io/en/latest/operations.html#spatial-wrapping

[tlogan2000]

Realigning all the longitudes of a netCDF is costly (0,360 to -180,180)

I think the ocgis wrapper would basically be doing this realignment of the dataset?

This is pretty easy to do in xarray but somewhat slow for really big datasets (The reason for why we chose to modify the bounds args instead of the data itself). Bascially we would need to first modify/adjust the xr.dataset lon coordinate values then call a xr.sortby() on the lon dim

One other question would be the expected output ...
I believe subset.polygon() in xclim currently wlll successfully subset the meridian crossing region but the resulting dataset still would have original 'all positive' lons so we end up with a region on the 'left' of a map figure and a region on the 'right' with nans inbetween

[Zeitsperre]

@tlogan2000 That's exactly how it works. The region is split and treated as multiple features. At the end of the day, since they're being used to build a masking layer, it doesn't really matter if there are multiple parts/features.

I'd like to see this type of functionality built into xarray directly since it's starting to feel like it's beyond the scope of xclim (calculating climate indices), but the projects to build out gdal/rasterio-like methods for xarray objects are still very new. Would it be worthwhile reopening this discussion and decide whether there's a project in particular we would want to back? (My recommendation goes to rioxarray or regionmask, but I would need to re-acquaint myself with recent changes)

[Zeitsperre]

Relevant issue concerning this issue for clisops.core.subset.subset_bbox → roocs/clisops#119