来自 netcdf 的 Xarray 数据数组返回大于输入的 numpy 网格数组

Question

我有一个 netcdf 文件，其中的浮点值表示经纬度上的叶绿素浓度。我试图在两组纬度/经度之间画一条线，并从线上的点返回所有叶绿素值。

我从几何学的角度接近它：对于点 (x1, y1) 和 (x2, y2)，找到直线的斜率和截距，并返回给定 y 值的所有 x 值线。一旦我有了所有的 x 和 y 值（经度和纬度），我希望将它们输入到 xarray 选择方法中以返回叶绿素浓度。

ds = '~/apr1.nc'
ds = xarray.open_dataset(ds, decode_times=False)

x1, y1 = [34.3282, 32.4791]
x2, y2 = [34.7, 32.21]

slope = (y2 - y1) / (x2 - x1)
intercept = y1 - (slope * x1)

line_lons = np.arange(x1, x2, step)
line_lats = [slope * x + intercept for x in lons]
values = ds.CHL.sel(lat=line_lats, lon=line_lons, method='nearest')
ds.values

>>> [0.0908799 , 0.06634101, 0.07615771, 0.16289435],
            [0.06787204, 0.07480557, 0.0655338 , 0.06064864],
            [0.06352911, 0.06586582, 0.06702182, 0.10024723],
            [0.0789495 , 0.07035938, 0.07455409, 0.08405576]]], dtype=float32)

line_lons
>>> array([34.3282, 34.4282, 34.5282, 34.6282])

我想创建一个以 x 轴为经度、以 y 轴为值的绘图。问题是 ds.values 命令返回一个形状为 (1, 4, 4) 的 numpy 数据数组，而经度只有 4。返回的数组中有更多的值。

plt.plot(line_lons, chlvalues.values)

知道为什么会这样以及如何为一个输入返回一个值吗？

谢谢。

Answer 1

我认为这是因为默认情况下您的输出是从方框中获取的，而不是沿着选定的样带。

我建议使用Numpy 和netCDF4 提供更复杂的解决方案，您首先使用随机坐标制作横断面，然后将这些随机坐标转换为输入文件中最接近的唯一坐标（唯一 = 以便沿样带只计算一次）。

之后，当你知道你的输出坐标时，你有两种可能性如何沿样带取出数据：

a) 你找到对应坐标的索引 b) 将原始数据插值到这些坐标（最近或双线性方法）

代码如下：

#!/usr/bin/env ipython
# --------------------------------------------------------------------------------------------------------------
import numpy as np
from netCDF4 import Dataset
# -----------------------------
# coordinates:
x1, y1 = [10., 55.]
x2, y2 = [20., 58.]
# --------------------------------
# ==============================================================================================================
# create some test data:
nx,ny = 100,100
dataout = np.random.random((ny,nx));
# -------------------------------
lonout=np.linspace(9.,30.,nx);
latout=np.linspace(54.,66.,ny);
# make data:
ncout=Dataset('test.nc','w','NETCDF3_CLASSIC');
ncout.createDimension('lon',nx);
ncout.createDimension('lat',ny);
ncout.createDimension('time',None);
ncout.createVariable('lon','float64',('lon'));ncout.variables['lon'][:]=lonout;
ncout.createVariable('lat','float64',('lat'));ncout.variables['lat'][:]=latout;
ncout.createVariable('var','float32',('lat','lon'));ncout.variables['var'][:]=dataout;
ncout.close()
#=================================================================================================================
# CUT THE  DATA FROM FILE:
# make some arbitrary line between start-end point, later let us convert it to indices:
coords=np.linspace(x1+1j*y1,x2+1j*y2,1000);
xo=np.real(coords);yo=np.imag(coords);
# ------------------------------------------------------
# get transect:
ncin = Dataset('test.nc');
lonin=ncin.variables['lon'][:];
latin=ncin.variables['lat'][:];
# ------------------------------------------------------
# get the transect indices:
rxo=np.array([np.squeeze(np.min(lonout[np.where(np.abs(lonout-val)==np.abs(lonout-val).min())])) for val in xo]);
ryo=np.array([np.squeeze(np.min(latout[np.where(np.abs(latout-val)==np.abs(latout-val).min())])) for val in yo]);
rcoords=np.unique(rxo+1j*ryo);
rxo=np.real(rcoords);ryo=np.imag(rcoords);
# ------------------------------------------------------
ixo=[int(np.squeeze(np.where(lonin==val))) for val in rxo];
jxo=[int(np.squeeze(np.where(latin==val))) for val in ryo];
# ------------------------------------------------------
# get var data along transect:
trans_data=np.array([ncin.variables['var'][jxo[ii],ixo[ii]] for ii in range(len(ixo))]);
# ------------------------------------------------------
ncin.close()
# ================================================================================================================
# Another solution using interpolation, when we already know the target coordinates (original coordinates along the transect):
from scipy.interpolate import griddata
ncin = Dataset('test.nc');
lonin=ncin.variables['lon'][:];
latin=ncin.variables['lat'][:];
varin=ncin.variables['var'][:];
ncin.close()
# ----------------------------------------------------------------------------------------------------------------
lonm,latm = np.meshgrid(lonin,latin);
trans_data_b=griddata((lonm.flatten(),latm.flatten()),varin.flatten(),(rxo,ryo),'nearest')