reproject icon indicating copy to clipboard operation
reproject copied to clipboard

Chess-like pattern when using rountrip_coords with GWCS, when using dask

Open astrofrog opened this issue 10 months ago • 2 comments

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:

chess

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.

astrofrog avatar Apr 09 '24 13:04 astrofrog

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:

chess

astrofrog avatar Apr 09 '24 13:04 astrofrog

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?

astrofrog avatar Apr 09 '24 13:04 astrofrog

Fixed now the file has been regenerated with https://github.com/astropy/reproject/pull/439

astrofrog avatar Jul 26 '24 10:07 astrofrog