xcube-dev / xcube

xcube is a Python package for generating and exploiting data cubes powered by xarray, dask, and zarr.
https://xcube.readthedocs.io/
MIT License
201 stars 20 forks source link

Resampling in space from EPGS:4326 (WGS84) to EPSG:32632 (UTM 32N) induces a shift #1073

Closed konstntokas closed 2 months ago

konstntokas commented 2 months ago

Describe the bug When using the function resampling_in_space() to project Sentinel-2 data retrieved with xcube-sh from EPGS:4326 (WGS84) to EPSG:32632 (UTM 32N), a shift has been noticed, when comparing to the reprojection method in rasterio.

To Reproduce Steps to reproduce the behavior with face data.

import matplotlib.pyplot as plt
import numpy as np
import pyproj
import xarray as xr

from xcube.core.resampling import resample_in_space
from xcube.core.gridmapping import GridMapping

# initialize dummy data
size = 1000
ds = xr.Dataset()
ds = ds.assign_coords(
    coords=dict(lon=np.linspace(9.6, 9.8, size), lat=np.linspace(47.6, 47.4, size))
)
data = (
    np.ones((size, size))
    * np.linspace(0, 10, size)
    * np.linspace(0, 10, size)[:, np.newaxis]
)
ds["B03"] = xr.DataArray(data, dims=("lat", "lon"))
print(ds)

# calculate transformed edge points
epsg_target = "EPSG:32632"
transformer = pyproj.Transformer.from_crs("EPSG:4326", epsg_target, always_xy=True)
psw = transformer.transform(ds.lon[0].values, ds.lat[0].values)
pse = transformer.transform(ds.lon[-1].values, ds.lat[0].values)
pne = transformer.transform(ds.lon[-1].values, ds.lat[-1].values)
pnw = transformer.transform(ds.lon[0].values, ds.lat[-1].values)
px = [psw[0], pse[0], pnw[0], pne[0]]
py = [psw[1], pse[1], pnw[1], pne[1]]

# perform reprojection via xcube
source_gm = GridMapping.from_dataset(ds)
temp_utm_target_gm = source_gm.transform(epsg_target)
utm_target_gm = temp_utm_target_gm.to_regular()
ds_reprojected = resample_in_space(ds, source_gm=source_gm, target_gm=utm_target_gm)
print(ds_reprojected)

fig, ax = plt.subplots(figsize=(7, 6))
ds_reprojected.B03.plot()
ax.plot(px, py, marker="x", linestyle="None", color="red")
plt.show()
konstntokas commented 2 months ago

Solution: When non-geographic target GridMapping is created in utm_target_gm = temp_utm_target_gm.to_regular(), the parameter is_lon_360 is set to true, because the values in x are above 180. Then 360 is subtracted from the x-axis to project the x-axis from [0, 360] -> [-180, 180], which induces a shift along the x-axis. This however should not be done for non-geographic target grids.