通过 3-D xarray.Dataset 绘制切片

Question

我想通过 3-D xarray.Dataset 绘制一个切片，如下所示：

这就是我开始的方式（为简单起见，使用教程air_temperature数据集）

# import packages
import xarray as xr
import numpy as np

# load dataset
ds = xr.tutorial.open_dataset("air_temperature")

# get slice values
tgt_lon = xr.DataArray(np.linspace(220, 280, num=15), dims="lon")
tgt_lat = xr.DataArray(np.linspace(30, 50, num=15), dims="lat")

# crop to region of interest - this works fine
da = ds.sel(lon=tgt_lon, 
            lat=tgt_lat, 
            method="nearest")

现在我们仍然有一个 3-D xarray.Dataset。对于二维绘图，我们希望将经度和纬度堆叠到相对于零点的距离数组（此处：角点 x0/y0，如图所示）。

# zero point: lat_min/lon_min
lon_orig = da.lon.min().values
lat_orig = da.lat.min().values

# stack longitude and latitude -- gives tuple for dist
sta = da.stack(dist=('lon', 'lat')) 

...然后我们可以通过循环计算相对于 lon_orig/lat_orig 的距离sta。对于实际绘图（使用contourf），我们会将数组展平为一维数组。这一切看起来有点乏味，我想知道我是否错过了这样做的明显方法？

最终，我们需要三个具有相同形状的一维数组来进行绘图：

• z（对于该数据集：时间）
• 距离（相对于 x0/y0，这里根据经度和纬度计算）
• 实际数据（此处：气温）

谢谢！

Answer 1

设置我的示例

首先，我将使用一个不同的例子，因为我想知道我应该得到什么结果。我的示例是一个 100x100x100 数组，其中包含一个半径为 50 的球。

我采用了这个答案中的代码来绘制它。

n = 100

def distance_from_center(x,y,z):
    V = np.stack([x-n/2,y-n/2,z-n/2])
    return np.linalg.norm(V, axis=0)

ball = np.fromfunction(distance_from_center, (n, n, n), dtype='float')
ball = (ball > 50).astype('int')

基本解决方案

只需 4 行代码就可以完成这样的事情。所以在我变得更喜欢之前我会先这样做。我的第一片是从一个角到另一个角。像这样：

我已经用绿色标记了要切片的位置。

t = np.arange(n)
slice = ball[t,t,:]
plt.matshow(slice.T, cmap='gray')
plt.gca().set_aspect(1/2**(1/2))

请注意，这里我必须将长宽比设置为二分之二的平方根，因为对角线方向上的 50 个像素的距离比普通方向上的 50 个像素的距离更长。图片看起来和预期的一样。两边都有空白，因为数组的角为其提供了空间，但我们仍然得到一个圆圈。

您可以像这样进行其他切割，但您总是在考虑如何获得整数坐标。

更复杂的解决方案

另一种方法是仅采用任何坐标并在未完美命中点的点之间进行插值。

import itertools

class Image_knn():
    def fit(self, image):
        self.image = image.astype('float')

    def predict(self, x, y, z):
        image = self.image
        weights_x = [(1-(x % 1)).reshape(-1), (x % 1).reshape(-1)]
        weights_y = [(1-(y % 1)).reshape(-1), (y % 1).reshape(-1)]
        weights_z = [(1-(z % 1)).reshape(-1), (z % 1).reshape(-1)]
        start_x = np.floor(x)
        start_y = np.floor(y)
        start_z = np.floor(z)

        return sum([image[np.clip(np.floor(start_x + x), 0, image.shape[0]-1).astype('int'),
                          np.clip(np.floor(start_y + y), 0, image.shape[1]-1).astype('int'),
                          np.clip(np.floor(start_z + z), 0, image.shape[1]-1).astype('int')]
                    *weights_x[x]*weights_y[y]*weights_z[z]
                    for x,y,z in itertools.product(range(2),range(2),range(2))])

    
image_model = Image_knn()
image_model.fit(ball)

fig, ax = plt.subplots(nrows=3,ncols=3)
ax = ax.reshape(-1)

我将给出一个使用不同切片的示例，其中一个方向是 z 方向，另一个方向是包含左下角和右下角变化点的方向。首先是指示切割位置的 3D 图：

现在我计算切割的坐标并绘制它们。这次我可以保持纵横比，因为我在较长的方向上使用了更多的点。请注意，我将第二个点的坐标放在绘图的标题中。

fig, ax = plt.subplots(nrows=3,ncols=3)
ax = ax.reshape(-1)

for i,x in enumerate(np.linspace(0,100,9)):
    p = np.array([0,100])
    q = np.array([100,100-x])
    distance = np.round(np.linalg.norm(q-p)).astype('int')
    t = np.linspace(0,1,distance)
    xy = t.reshape(-1,1)*q+(1-t).reshape(-1,1)*p
    r,s = np.meshgrid(np.arange(100), np.arange(distance))
    x = xy[s][:,:,0].reshape(-1)
    y = xy[s][:,:,1].reshape(-1)
    z = r.reshape(-1)

    out = image_model.predict(x,y,z)
    ax[i].imshow(out.reshape(distance, 100).T, cmap='gray',vmin=0,vmax=1)
    ax[i].set_title(tuple(q.astype('int')))