Xarray滚动和Cartopy add_cyclic_point函数冲突

Question

最近，我试着把180W-20E的太平洋和大西洋画成一张简单的图。为了跨越0度线，我使用了滚动方法，如下所示：

import matplotlib.pyplot   as plt
import cartopy.crs         as ccrs
import cartopy.feature     as cf
import cartopy             as cartopy
import numpy               as np
import xarray              as xr

from cartopy.mpl.ticker    import LongitudeFormatter, LatitudeFormatter
from cartopy.util          import add_cyclic_point


dset  =  xr.open_dataset('sst.mnmean.nc')
var   =  dset['sst'][:,:,:]
lat   =  dset['lat'][:]
lon   =  dset['lon'][:]

使用滚动方法：

lonr = lon.roll(lon=300)

这将为我提供以下经度数组：

[ 120.  122.  124.  126.  128.  130.  132.  134.  136.  138.  140.  142.
  144.  146.  148.  150.  152.  154.  156.  158.  160.  162.  164.  166.
  168.  170.  172.  174.  176.  178.  180.  182.  184.  186.  188.  190.
  192.  194.  196.  198.  200.  202.  204.  206.  208.  210.  212.  214.
  216.  218.  220.  222.  224.  226.  228.  230.  232.  234.  236.  238.
  240.  242.  244.  246.  248.  250.  252.  254.  256.  258.  260.  262.
  264.  266.  268.  270.  272.  274.  276.  278.  280.  282.  284.  286.
  288.  290.  292.  294.  296.  298.  300.  302.  304.  306.  308.  310.
  312.  314.  316.  318.  320.  322.  324.  326.  328.  330.  332.  334.
  336.  338.  340.  342.  344.  346.  348.  350.  352.  354.  356.  358.
    0.    2.    4.    6.    8.   10.   12.   14.   16.   18.   20.   22.
   24.   26.   28.   30.   32.   34.   36.   38.   40.   42.   44.   46.
   48.   50.   52.   54.   56.   58.   60.   62.   64.   66.   68.   70.
   72.   74.   76.   78.   80.   82.   84.   86.   88.   90.   92.   94.
   96.   98.  100.  102.  104.  106.  108.  110.  112.  114.  116.  118.]

然后我继续使用'sel‘方法提取我想要的区域：

lonr = lonr.sel(lon=slice(180,12))
lonr = lonr.roll(lon=104)

正如here所讨论的，最后一条线是为了避免在绘图上出现奇怪的水平线。也就是说，我的分片经度数组现在看起来像这样：

[   0.    2.    4.    6.    8.   10.   12.  180.  182.  184.  186.  188.
  190.  192.  194.  196.  198.  200.  202.  204.  206.  208.  210.  212.
  214.  216.  218.  220.  222.  224.  226.  228.  230.  232.  234.  236.
  238.  240.  242.  244.  246.  248.  250.  252.  254.  256.  258.  260.
  262.  264.  266.  268.  270.  272.  274.  276.  278.  280.  282.  284.
  286.  288.  290.  292.  294.  296.  298.  300.  302.  304.  306.  308.
  310.  312.  314.  316.  318.  320.  322.  324.  326.  328.  330.  332.
  334.  336.  338.  340.  342.  344.  346.  348.  350.  352.  354.  356.
  358.]

现在的问题是:这解决了区域选择问题，但它在绘图时给出了围绕0经度值的奇怪的白线，就像你可以看到的here一样。为了解决这个问题，我们会考虑使用"add_cyclic_point“方法，但这会产生以下错误：

var  =  var.sel(lat=slice(-30,30),
                lon=lonr)

lon_idx = var.dims.index('lon')
var_c, lon_c = add_cyclic_point(var.values, 
                                coord=lonr, axis=lon_idx)

(...)in add_cyclic_point
    raise ValueError('The coordinate must be equally spaced.')
ValueError: The coordinate must be equally spaced.

正如它所说的，我不能将点添加到数组中，因为它的间距不相等。从方法手册看，我不能使用它的原因看起来很清楚，但我当然不喜欢在我的图(an example of what I did can be seen here)上到处都是那条奇怪的白线。

很抱歉问了这么长的问题，但是有人知道如何解决这个问题吗？任何帮助都将不胜感激。

Answer 1

考虑到我们无法访问您的代码来重现问题，我将发明一些数字来说明以您所拥有的方式滚动的问题：

In [1]: import numpy as np

In [2]: lons = np.arange(0, 360, 10)

In [3]: lons
Out[3]: 
array([  0,  10,  20,  30,  40,  50,  60,  70,  80,  90, 100, 110, 120,
       130, 140, 150, 160, 170, 180, 190, 200, 210, 220, 230, 240, 250,
       260, 270, 280, 290, 300, 310, 320, 330, 340, 350])

现在，如果我们像你所做的那样滚动：

In [4]: lons = np.roll(lons, -list(lons).index(100))

In [5]: lons
Out[5]: 
array([100, 110, 120, 130, 140, 150, 160, 170, 180, 190, 200, 210, 220,
       230, 240, 250, 260, 270, 280, 290, 300, 310, 320, 330, 340, 350,
         0,  10,  20,  30,  40,  50,  60,  70,  80,  90])

这些数字的问题是，如果你认为它们是球形的( 350离0只有10度)，那么它们就很好，但是如果你认为它们是笛卡尔的，那么有一个点，你的两个数据点之间的距离是350度。绘制数据时，如果您告诉cartopy您的数据在PlateCarree投影坐标系中，这正是正在发生的事情。

在这种情况下，有许多方法可以固定您的数字，以便您的一些点之间不会有如此大的距离。一种方法：

In [6]: np.diff(lons) < -180
Out[6]: 
array([False, False, False, False, False, False, False, False, False,
       False, False, False, False, False, False, False, False, False,
       False, False, False, False, False, False, False,  True, False,
       False, False, False, False, False, False, False, False], dtype=bool)

In [7]: np.cumsum(np.diff(lons) < -180)
Out[7]: 
array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
       0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1])

In [8]: new_lons = lons.copy()

In [8]: new_lons[1:] += np.cumsum(np.diff(lons) < -180) * 360

In [9]: new_lons
Out[9]: 
array([100, 110, 120, 130, 140, 150, 160, 170, 180, 190, 200, 210, 220,
       230, 240, 250, 260, 270, 280, 290, 300, 310, 320, 330, 340, 350,
       360, 370, 380, 390, 400, 410, 420, 430, 440, 450])

您现在回到了投影空间中等距坐标的领域，并且应该更幸运地使用像add_cyclic_point这样的东西。