Объединение наборов данных xarray с одинаковым экстентом, но разным пространственным разрешением

Question

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

Честно говоря, я не уверен, что это лучший способ сравнить эти два значения при поиске влияния осадка на SIF, поэтому я открыт для идей различных подходов. Что касается объединения данных в настоящее время, я использую xr.combine_by_coords, но это дает мне ошибку, которую я описал ниже. Я мог бы также сделать это, преобразовав netcdfs в геотифов, а затем используя растерио, чтобы деформировать их, но это, похоже, неэффективный способ сделать это сравнение. Вот что у меня пока есть:

import netCDF4
import numpy as np
import dask
import xarray as xr

rainy_bbox = np.array([
    [-69.29519955115512,-13.861261028444734],
    [-69.29519955115512,-12.384786628185896],
    [-71.19583431678012,-12.384786628185896],
    [-71.19583431678012,-13.861261028444734]])
max_lon_lat = np.max(rainy_bbox, axis=0)
min_lon_lat = np.min(rainy_bbox, axis=0)

# this dataset is available here: ftp://fluo.gps.caltech.edu/data/tropomi/gridded/
sif = xr.open_dataset('../data/TROPO_SIF_03-2018.nc')

# the dataset is global so subset to my study area in the Amazon
rainy_sif_xds = sif.sel(lon=slice(min_lon_lat[0], max_lon_lat[0]), lat=slice(min_lon_lat[1], max_lon_lat[1]))  

# this data can all be downloaded from NASA Goddard here either manually or with wget but you'll need an account on https://disc.gsfc.nasa.gov/: https://pastebin.com/viZckVdn
imerg_xds = xr.open_mfdataset('../data/3B-DAY.MS.MRG.3IMERG.201803*.nc4')

# spatial subset
rainy_imerg_xds = imerg_xds.sel(lon=slice(min_lon_lat[0], max_lon_lat[0]), lat=slice(min_lon_lat[1], max_lon_lat[1]))  

# I'm not sure the best way to combine these datasets but am trying this
combo_xds = xr.combine_by_coords([rainy_imerg_xds, rainy_xds])

В настоящее время я получаю, казалось бы, бесполезную RecursionError: maximum recursion depth exceeded in comparison в этой последней строке. Когда я добавляю аргумент join='left', тогда данные из набора данных rainy_imerg_xds находятся в combo_xds, а когда я делаю join='right', данные rainy_xds присутствуют, а если я делаю join='inner', то никаких данных нет. Я предположил, что с этой функцией была некоторая внутренняя интерполяция, но, похоже, нет.

Answer 1

Эта документация от xarray довольно просто описывает решение этой проблемы. xarray позволяет интерполировать в нескольких измерениях и задавать измерения x и y другого набора данных в качестве выходных измерений. Таким образом, в этом случае это делается с

# interpolation based on http://xarray.pydata.org/en/stable/interpolation.html
# interpolation can't be done across the chunked dimension so we have to load it all into memory
rainy_sif_xds.load()

#interpolate into the higher resolution grid from IMERG
interp_rainy_sif_xds = rainy_sif_xds.interp(lat=rainy_imerg_xds["lat"], lon=rainy_imerg_xds["lon"])

# visualize the output
rainy_sif_xds.dcSIF.mean(dim='time').hvplot.quadmesh('lon', 'lat', cmap='jet', geo=True, rasterize=True, dynamic=False, width=450).relabel('Initial') +\
interp_rainy_sif_xds.dcSIF.mean(dim='time').hvplot.quadmesh('lon', 'lat', cmap='jet', geo=True, rasterize=True, dynamic=False, width=450).relabel('Interpolated')

# now that our coordinates match, in order to actually merge we need to convert the default CFTimeIndex to datetime to merge dataset with SIF data because the IMERG rainfall dataset was CFTime and the SIF was datetime
rainy_imerg_xds['time'] = rainy_imerg_xds.indexes['time'].to_datetimeindex()

# now the merge can easily be done with
merged_xds = xr.combine_by_coords([rainy_imerg_xds, interp_rainy_sif_xds], coords=['lat', 'lon', 'time'], join="inner")

# now visualize the two datasets together // multiply SIF by 30 because values are so ow
merged_xds.HQprecipitation.rolling(time=7, center=True).sum().mean(dim=('lat', 'lon')).hvplot().relabel('Precip') * \
(merged_xds.dcSIF.mean(dim=('lat', 'lon'))*30).hvplot().relabel('SIF')