Xarray Roll 和 Cartopy add_cyclic_point 函数冲突

Question

最近我尝试制作太平洋和大西洋的简单图，从 180W 到 20E。为了越过 0 度线，我使用了 roll 方法如下：

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 之类的东西有更多的运气。