simonsobs / Pixell.jl

next-generation sky map manipulation on rectangular pixels
MIT License
3 stars 1 forks source link

memory copy issue when using sub on maps with more than 2 dims #42

Closed xzackli closed 2 years ago

xzackli commented 2 years ago

I think the sub function is affecting some memory references in the C struct, and does bad things to the WCSTransform.

This manifests in a strange way: the copy method is broken on such WCS.

Here's the MWE @guanyilun

using Pixell, WCS

shape, wcs = fullsky_geometry(WCSTransform, deg2rad(1))
shape = (shape..., 2)  # this is a map with shape (nx, ny, 2)

# read map to disk and then get it back
m = Enmap(randn(shape), wcs)
write_map("test.fits", m)
m2 = read_map("test.fits")  # add ; trim=false if PR #41 is merged
wcs = m2.wcs 

# check sky coordinates, then check sky coordinates of a *copy*
print(pix_to_world(wcs, [50., 50.]), " ", pix_to_world(copy(wcs), [50., 50.]))
[131.0, -41.0] [5.45352918278e-312, 0.0]

This is related to https://github.com/JuliaAstro/WCS.jl/issues/43

guanyilun commented 2 years ago

ok, I think I've figured out what is going on. When the image is loaded from disk, sub is called to trim the wcs to leave only the first two axes. I did it in a sloppy way by simply changing the naxis field. This works for all 1D fields such as crval, cdelt, etc., but it doesn't change naxis and the linear transformation matrix stored in the linprm struct in wcs.lin. For example, after the map has been read from disk (m2), wcs.naxis=2 but wcs.lin.naxis=3. This mismatch causes memory overflow when one tries to make a copy of it, resulting in the above error.

I fixed it in #47 by rewriting sub using wcslib.