Xarray Data Array из netcdf возвращает массив массивов, превышающий ввод

Question

У меня есть файл netcdf со значениями с плавающей запятой, представляющими концентрацию хлорофилла в широтах и ​​долготах. Я пытаюсь нарисовать линию между двумя наборами латов и долг и вернуть все значения хлорофилла из точек на линии.

Я подхожу к нему с геометрической точки зрения: для точек (x1, y1) и (x2, y2) найдите наклон и точку пересечения линии и верните все значения x для заданных значений y на линия. Как только у меня будут все значения 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), в то время как долготы равны 4. В возвращаемом массиве значений намного больше.

plt.plot(line_lons, chlvalues.values)

Есть идеи, почему это так и как я могу вернуть одно значение для одного входа?

Спасибо.

Answer 1

Я предполагаю, что это потому, что по умолчанию ваш вывод был взят из коробки, а не вдоль выбранного разреза.

Я предлагаю более сложное решение с Numpy и netCDF4, где вы сначала создаете трансект со случайными координатами, а затем превращаете эти случайные координаты в ближайшие ближайшие уникальные координаты из входного файла (уникальные =, чтобы каждая точка вдоль трансект включается только один раз).

Впоследствии, когда вы знаете свои выходные координаты, у вас есть 2 возможности вывести данные по разрезу:

а) вы найдете индексы соответствующих координат б) интерполировать исходные данные по этим координатам (ближайший или билинейный метод)

Вот код:

#!/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')