bug in `xcube.core.resampling.resample_in_space when upsampling between two different EPSGs

[konstntokas]

Describe the bug
When loading Sentinel-2 data with EPSG:3035 and CCI and C3S land cover classification map provided from DeepESDL public AWS S3 bucket, and trying to resample the land cover classification map given in EPSG:4326 to the the fine Sentinel-2 data in EPSG:3035 returns a land cover map filled with zeros.

To Reproduce
Here is an example which produces the bug:

import os

os.environ["NUMBA_DISABLE_JIT"] = "1"

from xcube.core.store import new_data_store
from xcube.core.resampling import resample_in_space
from xcube.core.gridmapping import GridMapping
import matplotlib.pyplot as plt
import pyproj
import datetime


# create sentinelhub datastore
store = new_data_store("sentinelhub")

# create s3 data store
s3_store = new_data_store(
    "s3",
    root="xxxx",
    storage_options=dict(
        anon=False,
        key="xxxx",
        secret="xxx",
    ),
)

# selection parameters
data_id = "S2L2A"
bbox = [9.6, 47.4, 9.8, 47.6]
tile_size = [1024, 1024]
time_range = ("2020-06-01", "2020-07-01")
variable_names = ["B03", "B04"]

# get data in EPSG:3035
transformer = pyproj.Transformer.from_crs("EPSG:4326", "EPSG:3035", always_xy=True)
pmin = transformer.transform(bbox[0], bbox[1])
pmax = transformer.transform(bbox[2], bbox[3])
bbox_3035 = [pmin[0], pmin[1], pmax[0], pmax[1]]
spatial_res_3035 = 20
ds_3035 = store.open_data(
    data_id,
    variable_names=variable_names,
    bbox=bbox_3035,
    spatial_res=spatial_res_3035,
    time_range=time_range,
    tile_size=tile_size,
    crs="EPSG:3035",
)
print(ds_3035)

# get land cover from S3
land_cover_levels = s3_store.open_data("LC-1x2160x2160-1.0.0.levels")
land_cover = land_cover_levels.get_dataset(0)
buffer = 0.3
land_cover_2019 = land_cover.sel(
    lat=slice(bbox[3] + buffer, bbox[1] - buffer),
    lon=slice(bbox[0] - buffer, bbox[2] + buffer),
    time=datetime.datetime(2019, 1, 1),
)

# resample the data
land_cover_2019_resamp = resample_in_space(land_cover_2019, ref_ds=ds_3035.isel(time=4))
print(land_cover_2019_resamp)
land_cover_2019_resamp.lccs_class.plot()
plt.show()

Expected behavior
The function resample_in_space should change the projection of the land cover classification map and resample it to the grid of the Sentinel-2 data.

Additional context
I started debugging it, with partial success:
The function rectify_dataset line 160 returns a ndarray filled with nans since source_gm and target_gm are still in the two different EPSGs and cannot be compared with each other. Before that, the source_gm is created from a dataset source_ds_subset in line 141, where the dataset source_ds_subset has the old coordinates for EPSG:4326, and carries the new coordinates in source_ds_subset.transformed_x and source_ds_subset.transformed_y.