reproject
reproject copied to clipboard
Chess-like pattern when using rountrip_coords with GWCS, when using dask
I am seeing a very strange issue. The following example tries to reproject a dataset with a GWCS to a slightly different GWCS:
import os
import asdf
from copy import deepcopy
from astropy import units as u
from reproject import reproject_interp
import matplotlib.pyplot as plt
aia = asdf.open(
os.path.join(
"/Users/tom/Code/Astropy/reproject/reproject/tests/data/aia_171_level1.asdf"
)
)
input_data = (aia["data"][...], aia["wcs"])
wcs_out = deepcopy(aia["wcs"])
wcs_out.forward_transform.offset_0 = -60.3123 * u.pix
wcs_out.forward_transform.offset_1 = -61.9422 * u.pix
shape_out = aia["data"].shape
array1 = reproject_interp(
input_data,
wcs_out,
shape_out=shape_out,
return_footprint=False,
)
array2 = reproject_interp(
input_data,
wcs_out,
shape_out=shape_out,
return_footprint=False,
return_type="dask",
block_size=(200, 200),
)
array2 = array2.compute(scheduler="synchronous")
fig = plt.figure()
ax1 = fig.add_subplot(1, 2, 1)
ax1.imshow(array1, origin="lower")
ax2 = fig.add_subplot(1, 2, 2)
ax2.imshow(array2, origin="lower")
fig.savefig('chess.png')
The result is:
This goes away if using rountrip_coords=False
so is related to the round-tripping. However this seems to happen only when using dask, even if just in synchronous mode. Furthermore the chess-like pattern is unrelated to the chunk size. Even if I set the chunk size to be larger than the image, so that there is only one chunk, the chess-like pattern is present.
Ok so this is not related to dask per se but to the use, I think, of SlicedLowLevelWCS
:
import os
import asdf
from copy import deepcopy
from astropy import units as u
from reproject import reproject_interp
import matplotlib.pyplot as plt
from astropy.wcs.wcsapi import BaseHighLevelWCS, SlicedLowLevelWCS
from astropy.wcs.wcsapi.high_level_wcs_wrapper import HighLevelWCSWrapper
aia = asdf.open(
os.path.join(
"/Users/tom/Code/Astropy/reproject/reproject/tests/data/aia_171_level1.asdf"
)
)
input_data = (aia["data"][...], aia["wcs"])
wcs_out = deepcopy(aia["wcs"])
wcs_out.forward_transform.offset_0 = -60.3123 * u.pix
wcs_out.forward_transform.offset_1 = -61.9422 * u.pix
shape_out = aia["data"].shape
slices = (slice(0, 128), slice(0, 128))
if isinstance(wcs_out, BaseHighLevelWCS):
low_level_wcs = SlicedLowLevelWCS(wcs_out.low_level_wcs, slices=slices)
else:
low_level_wcs = SlicedLowLevelWCS(wcs_out, slices=slices)
wcs_out_sub = HighLevelWCSWrapper(low_level_wcs)
array = reproject_interp(
input_data,
wcs_out_sub,
shape_out=shape_out,
return_footprint=False,
)
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.imshow(array, origin="lower")
fig.savefig("chess.png")
gives:
This might be because I'm not supposed to change the forward transform like this, as the backward transform would also need to be changed?
Fixed now the file has been regenerated with https://github.com/astropy/reproject/pull/439