Problem with dateline (longitude=0) for higher resolution global grids

[trondkr]

When I use the grid_global method for generating a 1x1 degree grid where I interpolate CMIP6 global climate models I keep experiencing a problem with data at longitude=0. This is very visible in the resulting maps. If I use a global grid with resolution 4x4 the problem goes away. I have tried a multitude of approaches to deal with the problem and I thought adding cartopy add_cyclic_point would surely fix it but no. So I am thinking that I do something wrong using the regridder which produces erroneous data at longitude=0.

I have a rather long notebook where I use CMIP6 data for calculating various properties and I plot the global results. The notebook is found here.

I define my global grid using : ds_out = xe.util.grid_global(1, 1)

And my regridder as:

def regrid_variable(varname, ds_in, ds_out, transpose=True):
    regridder = xe.Regridder(ds_in, ds_out, 'bilinear', reuse_weights=False, ignore_degenerate=True)
    regridder._grid_in = None
    regridder._grid_out = None
    if transpose:
        return regridder(ds_in[varname].T)
    else:
        return regridder(ds_in[varname])

To see the problem it is probably best to look at the notebook that has the images, but I am hoping you would see an obvious error in my definition of the global grid or the regridder. 
Thanks for your help. Cheers, Trond

[ahuang11]

Maybe applying this may fix it: https://scitools.org.uk/cartopy/docs/latest/cartopy/util/util.html

applicable example https://gist.github.com/darothen/c7560d8d19ffca90024c1f2df4927599

[trondkr]

@ahuang11 Thanks for your response. I already do use the add_cyclic_point to my code and problem is still there. I think I do something wrong with the interpolation.

indata_cyclic, lon_cyclic = add_cyclic_point(indata, coord=lon)

[ahuang11]

Did you try adding it before the interpolation?

On Tue, Jun 2, 2020, 3:26 PM Trond Kristiansen notifications@github.com wrote:

@ahuang11 https://github.com/ahuang11 Thanks for your response. I already do use the add_cyclic_point to my code and problem is still there. I think I do something wrong with the interpolation.

indata_cyclic, lon_cyclic = add_cyclic_point(indata, coord=lon)

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/JiaweiZhuang/xESMF/issues/101#issuecomment-637787229, or unsubscribe https://github.com/notifications/unsubscribe-auth/ADU7FFTJB6ATEF75XHR22ZDRUVOBBANCNFSM4NQ7QSJA .

[ahuang11]

Actually you probably want to set periodic=True in your regrid keyword args.

On Tue, Jun 2, 2020, 3:46 PM Andrew H huang.andrew12@gmail.com wrote:

Did you try adding it before the interpolation?

On Tue, Jun 2, 2020, 3:26 PM Trond Kristiansen notifications@github.com wrote:

@ahuang11 https://github.com/ahuang11 Thanks for your response. I already do use the add_cyclic_point to my code and problem is still there. I think I do something wrong with the interpolation.

indata_cyclic, lon_cyclic = add_cyclic_point(indata, coord=lon)

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/JiaweiZhuang/xESMF/issues/101#issuecomment-637787229, or unsubscribe https://github.com/notifications/unsubscribe-auth/ADU7FFTJB6ATEF75XHR22ZDRUVOBBANCNFSM4NQ7QSJA .

[trondkr]

@ahuang11 I tried that too without success: regridder = xe.Regridder(ds_in, ds_out, 'bilinear', reuse_weights=True, periodic=True, ignore_degenerate=False)

However, I did not try to add cyclic point to a dataset prior to interpolation. Is there an easy way to do that, which takes care of coordinates and variables in the dataset simultaneously?

[trondkr]

@ahuang11 So it seems my problem is caused by my input grid having larger steps between each grid point than I thought.

(y) float64 45.0 46.25 47.5 48.75 ... 86.25 87.5 88.75 90.0 (x) float64 0.0 1.875 3.75 5.625 ... 352.5 354.4 356.2 358.1

The longitude ends at 358.1. If I use a global grid with 2x2 degrees resolution I do not experience the line at longitude 359/0, but with a step size of 1 longitude, I do. Is there a good way to fix this without having to increase my grid to 2x2 instead of 1x1?

[trondkr]

I solved my problem by doing the interpolation twice (not the most elegant solution), first using a coarser resolution regridder to get global coverage and then downscaling the solution. Suggestions for improvements are very welcome. But it works...

ds_out_amon = xe.util.grid_2d(-180,180,2,southern_limit_latitude,90,2) 
ds_out = xe.util.grid_2d(-180,180,1,southern_limit_latitude,90,1) #grid_global(1, 1)

dr_out_uas_amon=regrid_variable("uas",ds_uas,ds_out_amon,transpose=True).to_dataset()
dr_out_uas=regrid_variable("uas",dr_out_uas_amon,ds_out,transpose=False)`

[brews]

For what it's worth, I have a simple example of using cartopy.util.add_cyclic_point to add a wrap-around or cyclic pixel to xarray.Dataset here: https://gist.github.com/brews/95e167fb7df997864304f07ee321f4fd.

This removes the longitude band artifact if you add the wrap-around pixels before passing the data fields to instantiate xesmf.Regridder.