Fixes issue #30 and is related to PR #28

[razman786]

Fixes issue https://github.com/ska-sa/tigger-lsm/issues/30 and is related to PR https://github.com/ska-sa/tigger-lsm/pull/28. Provides a self-contained Tigger v1.6.0 workaround.

[ratt-priv-ci]

Can one of the admins verify this patch?

[bennahugo]

No what I mean is we don't have a guarantee that axis 0 and axis 1 is the ra and dec axes. All that needs to be changed here is to use the axis id as read off the CTYPE

On Tue, May 11, 2021 at 2:37 AM Raz @.***> wrote:

@.**** commented on this pull request.

In Tigger/Coordinates.py https://github.com/ska-sa/tigger-lsm/pull/31#discussion_r629769841:

• use WCS pixel information to build refsky

• self.wcs = WCS(header)
• wcs_pixel = self.wcs.pixel_shape
• refsky_pixels = []
• for ipixel in range(naxis):
• if ipixel < 2:
• refsky_pixels.append(int(wcs_pixel[ipixel] / 2))
• else:
• refsky_pixels.append(wcs_pixel[ipixel])
• self.refsky = self.wcs.wcs_pix2world([refsky_pixels], 0)[0, :]

• set centre x/y pixels

• self.xpix0, self.ypix0 = self.refpix[self.ra_axis], self.refpix[self.dec_axis]

• set x/y scales

• pixscale = (header["CDELT1"], header["CDELT2"])

My first objective was to simply take the old tigger-lsm and use it with astropy instead of astLib, as a self contained workaround for the GUI that matches the old tigger-lsm (which it does). However, if this is unsafe, then its probably best to use something along the lines of:

from astropy.wcs import utilspixel_scales = utils.proj_plane_pixel_scales(self.wcs.celestial)self.xscale = -pixel_scales[0] DEGself.yscale = pixel_scales[1] DEG

The above outputs the correct scale for the Tigger GUI.

— You are receiving this because your review was requested. Reply to this email directly, view it on GitHub https://github.com/ska-sa/tigger-lsm/pull/31#discussion_r629769841, or unsubscribe https://github.com/notifications/unsubscribe-auth/AB4RE6WPCIF6VCXJOAPVEKLTNB35LANCNFSM44RL3BCA .

--

Benjamin Hugo

PhD. student, Centre for Radio Astronomy Techniques and Technologies Department of Physics and Electronics Rhodes University

Junior software developer Radio Astronomy Research Group South African Radio Astronomy Observatory Black River Business Park Observatory Cape Town

[razman786]

@bennahugo and @o-smirnov - FITSWCS has been re-implemented using Astropy methods, it also no longer inherits from FITSWCSpix. The new lm() and radec() methods should satisfy the related discussions for this PR.

This change came about from previous discussions with @bennahugo and that I discovered a test file, namely model/2015/combined-4M5S.fits, has NAXIS as 3 and WCS AXES as 4, pix2world then fails expecting N x 4.

Number of WCS axes: 4
CTYPE : 'RA---SIN'  'DEC--SIN'  'STOKES'  'FREQ'  
CRVAL : 85.65057465185  49.85200932178  1.0  1776500000.0  
CRPIX : 2049.0  2049.0  1.0  1.0  
PC1_1 PC1_2 PC1_3 PC1_4  : 1.0  0.0  0.0  0.0  
PC2_1 PC2_2 PC2_3 PC2_4  : 0.0  1.0  0.0  0.0  
PC3_1 PC3_2 PC3_3 PC3_4  : 0.0  0.0  1.0  0.0  
PC4_1 PC4_2 PC4_3 PC4_4  : 0.0  0.0  0.0  1.0  
CDELT : -0.0006944444444444  0.0006944444444444  1.0  384000000.0  
NAXIS : 4096  4096  4
WCS None
naxis 3

Using Astropy WCS methods and not sub-classing FITSWCSpix avoids the error and a reliance on the naxis. This error occurs on the old tigger-lsm and the current upstream version.

The offset() method has also been changed back to match SinWCS, as it seems to be working fine. It might be an idea to update the offset() method to also use Astropy calls - I did have a quick hack with it.

Re-implementations have been tested against the old version of tigger-lsm and the current upstream version, with Python 3 only.

This PR now has a standalone FITSWCS that works well with the new tigger. Comments and testing welcomed \o/

[bennahugo]

ok to test

[bennahugo]

Hi @razman786 if the naxis is not set to 4 (and does not conform to the shape of the stored datacube therefore) and the ctype4 is set to FREQ as it appears to be it indicates a non-conformant fits file to me. But either way it appears the WCS library applies a fix in this case

[bennahugo]

Alright this has passed however the ra dec coordinates tigger is printing is flipped as opposed to the shipped tigger lsm version. I'm still not convinced about either "SIN" implementations, but I will raise that as another issue as this is a preexisting issue

[bennahugo]

ie. your ra increases in the wrong direction @razman786 when running tigger 1.4

[bennahugo]

The following patch fixes this

From 928ab238970c5a1ee82c29a11bf534f34723adc0 Mon Sep 17 00:00:00 2001
From: bennahugo <bennahugo@gmail.com>
Date: Thu, 13 May 2021 13:19:53 +0200
Subject: [PATCH] Fix RA direction

---
 Tigger/Coordinates.py | 18 +++++++++---------
 1 file changed, 9 insertions(+), 9 deletions(-)

diff --git a/Tigger/Coordinates.py b/Tigger/Coordinates.py
index 5d8a33d..6f35dc1 100644
--- a/Tigger/Coordinates.py
+++ b/Tigger/Coordinates.py
@@ -252,7 +252,7 @@ class Projection(object):
                 refpix1[self.ra_axis] += 1
                 refpix1[self.dec_axis] += 1
                 delta = self.wcs.wcs_pix2world([refpix1], 0)[0] - self.refsky
-                self.xscale = delta[self.ra_axis] * DEG
+                self.xscale = -delta[self.ra_axis] * DEG
                 self.yscale = delta[self.dec_axis] * DEG
                 has_projection = True
             except Exception as exc:
@@ -265,7 +265,7 @@ class Projection(object):

         def lm(self, ra, dec):
             if not self.has_projection():
-                return numpy.sin(ra) / self.xscale, numpy.sin(dec) / self.yscale
+                return numpy.sin(ra) / -self.xscale, numpy.sin(dec) / self.yscale
             if numpy.isscalar(ra) and numpy.isscalar(dec):
                 if ra - self.ra0 > math.pi:
                     ra -= 2 * math.pi
@@ -295,7 +295,7 @@ class Projection(object):

         def radec(self, l, m):
             if not self.has_projection():
-                return numpy.arcsin(l * self.xscale), numpy.arcsin(m * self.yscale)
+                return numpy.arcsin(l * -self.xscale), numpy.arcsin(m * self.yscale)
             if numpy.isscalar(l) and numpy.isscalar(m):
                 pixvec = np.array(self.refpix).copy()
                 pixvec[self.ra_axis] = l
@@ -320,7 +320,7 @@ class Projection(object):

         def offset(self, dra, ddec):
             """ dra and ddec must be in radians """
-            return self.xpix0 + dra / self.xscale, self.ypix0 + ddec / self.xscale
+            return self.xpix0 + dra / -self.xscale, self.ypix0 + ddec / self.xscale
             # TODO - investigate; old code has 'self.xpix0 - dra...', new code is 'self.xpix0 + dra...'?
             # return self.xpix0 - dra / self.xscale, self.ypix0 + ddec / self.xscale

@@ -378,7 +378,7 @@ class Projection(object):

             # set x/y scales
             pix_scales = self.wcs.wcs.cdelt
-            self.xscale = pix_scales[self.ra_axis] * DEG
+            self.xscale = -pix_scales[self.ra_axis] * DEG
             self.yscale = pix_scales[self.dec_axis] * DEG

             # set l0, m0
@@ -396,12 +396,12 @@ class Projection(object):
                 l, m = -0.0, 0.0
             else:
                 l, m = coord_pixels[self.ra_axis], coord_pixels[self.dec_axis]
-            l = (l - self._l0) * self.xscale
+            l = (l - self._l0) * -self.xscale
             m = (m - self._m0) * self.yscale
             return l, m

         def radec(self, l, m):
-            x = self.xpix0 - l / self.xscale
+            x = self.xpix0 + l / -self.xscale
             y = self.ypix0 + m / self.yscale
             coord = utils.pixel_to_skycoord(xp=x, yp=y, wcs=self.wcs, origin=0, mode='all')
             ra = coord.ra.value
@@ -468,10 +468,10 @@ class Projection(object):

         def lm(self, ra, dec):
             l, m = Projection.FITSWCSpix.lm(self, ra, dec)
-            return sin((l - self._l0) * self.xscale), sin((m - self._m0)*self.yscale)
+            return sin((l - self._l0) * -self.xscale), sin((m - self._m0)*self.yscale)

         def radec(self, l, m):
-            return Projection.FITSWCSpix.radec(self, arcsin(l / self.xscale + self._l0), arcsin(m / self.yscale + self._m0))
+            return Projection.FITSWCSpix.radec(self, arcsin(l / -self.xscale + self._l0), arcsin(m / self.yscale + self._m0))

         def offset(self, dra, ddec):
             return sin(dra), sin(ddec)
-- 
2.17.1

[bennahugo]

I've tested this with tigger 1.6 and tigger 1.4 (py2.7)