Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

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

Closed
konstntokas opened this issue Sep 20, 2024 · 1 comment · Fixed by #1074
Closed

Comments

@konstntokas
Copy link
Contributor

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
Copy link
Contributor Author

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.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging a pull request may close this issue.

1 participant