Вычислить отклонение или аномалию значения между двумя массивами с разным размером географической сетки

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

Мои два файла netcdf имеют следующую структуру:

Наблюдения (инструментальные измерения):

Габаритные размеры:

lat: 1321, lon: 1321

Переменные данных:

Longitude (lon) float64 8.413 8.411 8.409 ... 4.904 4.905
Latitude (lat) float64 47.4 47.4 47.41 ... 52.37 52.37
obs_data (lat, lon) float64 ...

Данные модели:

Габаритные размеры:

latitude: 140, level: 1, longitude: 215, time: 24

Координаты:

longitude  (longitude)  float32    357.55 357.65 ... 18.85 18.95 
latitude   (latitude)   float32    55.95 55.85 55.75 ... 42.15 42.05    
level      (level)      float32    0.0
time       (time)    timedelta64[ns]    00:00:00 01:00:00 ... 23:00:00

Переменные данных:

model_data (time, level, latitude, longitude) float32 ...

Я пробовал всевозможные подходы, но каждый раз сталкиваюсь с какой-то ошибкой, для которой, похоже, нет решения, и мне приходилось пробовать другой подход. Самое близкое, что я получил, - это следование этому замечательному руководству, но и здесь я ударился о стену. Когда я пытаюсь найти ближайшие широты и долготы для двух наборов данных,

lat_idx = np.abs(model_lat - obs_lat).argmin() #subtract train lat from model lat
lon_idx = np.abs(model_lon - obs_lon).argmin() #subtract train lon from model lon

Я получаю следующую ошибку

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-437-9396b00ba22f> in <module>
     18 
     19 # Find the nearest latitude and longitude for the train data
---> 20 lat_idx = np.abs(model_lat - obs_lat).argmin()
     21 lon_idx = np.abs(model_lon - obs_lon).argmin()
     22 

~/opt/anaconda3/lib/python3.7/site-packages/numpy/ma/core.py in __sub__(self, other)
   4115         if self._delegate_binop(other):
   4116             return NotImplemented
-> 4117         return subtract(self, other)
   4118 
   4119     def __rsub__(self, other):

~/opt/anaconda3/lib/python3.7/site-packages/numpy/ma/core.py in __call__(self, a, b, *args, **kwargs)
   1024         with np.errstate():
   1025             np.seterr(divide='ignore', invalid='ignore')
-> 1026             result = self.f(da, db, *args, **kwargs)
   1027         # Get the mask for the result
   1028         (ma, mb) = (getmask(a), getmask(b))

ValueError: operands could not be broadcast together with shapes (140,) (1321,)

Разве нет способа просто посчитать:

anomaly = model_data[lat, lon] - obs_data[lat, lon]

?

Моя последняя надежда - xarray, но я действительно борюсь с их документацией, и я потратил дни на поиски пути вперед.

Кто-нибудь из вас нашел решение этой проблемы? Любые советы действительно приветствуются.

Редактировать:

По просьбе В. Айрата:

In: type(model_data)
Out: xarray.core.dataset.Dataset

obs_data - того же типа.

Если два значения obs_data попадают в одну и ту же model_data ячейку, obs_data следует вычесть из той же model_data ячейки.