Masking for observer

[christophwelling]

Fixes masking of entries below the horizon in the observer class Addresses issue #3

[codecov[bot]]

Codecov Report

Merging #4 (8b929e6) into master (0c98df2) will not change coverage. The diff coverage is 100.00%.

:exclamation: Current head 8b929e6 differs from pull request most recent head ef98725. Consider uploading reports for the commit ef98725 to get more accurate results

@@           Coverage Diff           @@
##           master       #4   +/-   ##
=======================================
  Coverage   90.09%   90.09%           
=======================================
  Files           8        8           
  Lines         333      333           
=======================================
  Hits          300      300           
  Misses         33       33           

Impacted Files	Coverage Δ
pygdsm/base_observer.py	96.55% <100.00%> (ø)

---

Continue to review full report at Codecov.

Legend - Click here to learn more Δ = absolute <relative> (impact), ø = not affected, ? = missing data Powered by Codecov. Last update 0c98df2...ef98725. Read the comment docs.

[telegraphic]

Hi @christophwelling, can you clarify what you think is wrong with the masking?

I generated 24-hr movies for observers at the North pole, on the equator, and at the South pole using the code in master:

(incidentally, placing an observer exactly at +/- 90 degrees causes issues!). This matches what I expect: one hemisphere (above the horizon) is visible, sources don't set if observing from the poles.

With the fix_observer_masking patch the observed skies look like this: which is not the sky view that I intended.

Perhaps I have misunderstood, or the underlying issue that brought this to your attention is elsewhere? (Coordinate transforms make my head hurt).

[christophwelling]

Then I am a bit confused about the coordinates. I thought theta was the elevation and phi the azimuth. Is it the other way around?

[telegraphic]

I am also a bit confused! Stepping through, we start in galactic coordinates, and applied the rotation

hrot = hp.Rotator(rot=[ra_deg, dec_deg], coord=['G', 'C'], inv=True)

Which transforms from galactic to equatorial coordinates, and applies a rotation to put RA, DEC at the center. The rotated coordinates are

g0, g1 = hrot(self._theta, self._phi)

I think that g0 and g1 are alt/az? I don't recall why the mask doesn't use g0 and g1, and instead uses _theta and _phi. But guessing that's where our confusion comes from?

[christophwelling]

self._theta and self._phi are taken from the healpy pix2ang function. If I understand that function correctly, this should mean theta is the co-latitude, so that one should be used to determine if something is below the horizon. I don't really understand exactly what the rotations do. I get the idea, but unfortunately healpy's documentation is kind of a mess, so it's hard for me to try to find any errors there.

[christophwelling]

So what's the status here? How do we proceed with this?

[telegraphic]

Well, I think the code output is correct, but that it is using the wrong coordinates?

Does this look better?

        # Apply rotation
        hrot = hp.Rotator(rot=[ra_deg, dec_deg], coord=['G', 'C'], inv=True)
        g0, g1 = hrot(self._theta, self._phi)
        pix0 = hp.ang2pix(self._n_side, g0, g1)
        sky_rotated = sky[pix0]

        # Generate a mask for below horizon
        mask1 = g0 + np.pi / 2 > 2 * np.pi
        mask2 = g0 < np.pi / 2
        mask = np.invert(np.logical_or(mask1, mask2))

(Will need to test it, but just showing that I think we should be using the rotated coords, not original?)

[telegraphic]

Update: that ^^ doesn't look right to me in a quick test.

I do think it would be good to understand what's going on with coordinates so I don't think we should close this issue! Might need to get some input from a healpix guru...

[christophwelling]

Do you know any healpy gurus? :-) Unfortunately, the healpy documentation is a mess and borderline unusable. But just so I get this right: theta and phi are supposed to be zenith angle (with upwards theta=0 and the horizontal theta=90deg) and the azimuth, right? Shouldn't the mask be based on the values of theta? It doesn't make sense to me to use phi then or g0 (since that seems to be in galactic coordinates)

[christophwelling]

Here is how I understand that part of the code: hrot transforms from horizontal to galactic coordinates (it has inv=True). Then theta and phi are zenith/azimuth coordinates in the local coordinate system. Using hrot they are transformed to the galactic coordinates g0 and g1. Then ang2pix gives the correct pixel index and the value of that pixel is returned. This would mean that theta is the zenith angle and is what should be used for the masks.

[christophwelling]

I played around with this a bit and it looks like these changes don't work correctly. On the master the debug skymaps look right for the locations I am looking at. Here is where I get the problem though: I want to get the radio temperature for a given zenith and azimuth angle. So I use this:

radio_sky = sky_observer.generate(frequency) i_pix = healpy.pixelfunc.ang2pix(sky_observer._n_side, zenith, azimuth, lonlat=False, nest=True) temp = radio_sky[i_pix]

I would expect temp to return a masked value for zenith>pi, but that is not what I get. Instead the masking looks pretty random. Am I doing something wrong here? How can I get the values for a specific direction?

[cdilullo]

This masking issue is still not solved, right? Because I'm interested in getting this solved so I can use PyGDSM to generate sky maps for doing astronomical-based calibration of the LWA. I think I may be on to something in terms of figuring out the HEALPix weirdness