Add elevation limit masking for sky coverage maps

[telegraphic]

Some gist code:

import pygsm
import ephem
import healpy as hp

sky = pygsm.GlobalSkyModel()
d = sky.generate(1420)
d = hp.ma(d)

## Parkes observatory
# * Lat: -33
# * Long: 148
# * Elevation limit: 30.5

th, ph = hp.pix2ang(512, np.arange(hp.nside2npix(512)))
hrot = hp.Rotator(coord=['G', 'C'], inv=True)
th, ph = hrot(th, ph)
pix = hp.ang2pix(512, th, ph)
dd = d[pix]

th0, ph0 = hp.pix2ang(512, np.arange(hp.nside2npix(512)))
dec_limit = np.deg2rad(-33 + 90 - 30.5) # lat + pi/2 - elev limit 
th_limit = th0 > np.pi/2 - dec_limit
dd.mask = np.invert(th_limit)
dec_limit

hp.mollview(np.log2(dd), coord=['C', 'G'], 
            cmap='gist_heat', title='Parkes observable sky')
hp.mollview(np.log2(dd), coord=['C'], 
            cmap='gist_heat', title='Parkes observable sky')

## Green Bank telescope
# * Lat: 38
# * Long: 79.5
# * Elevation limit: 5 degrees

d = sky.generate(1420)
d = hp.ma(d)

th, ph = hp.pix2ang(512, np.arange(hp.nside2npix(512)))
hrot = hp.Rotator(coord=['G', 'C'], inv=True)
th, ph = hrot(th, ph)
pix = hp.ang2pix(512, th, ph)
dd = d[pix]

th0, ph0 = hp.pix2ang(512, np.arange(hp.nside2npix(512)))
dec_limit = np.deg2rad(38 + 90 - 5 ) # lat + pi/2 - elev limit 
th_limit = th0 < dec_limit #- np.deg2rad(15)
dd.mask = np.invert(th_limit)
dec_limit

#hp.mollview(th_limit)
hp.mollview(np.log2(dd), coord=['C', 'G'], 
            cmap='inferno', title='Green Bank observable sky')
hp.mollview(np.log2(dd), coord=['C'], 
            cmap='inferno', title='Green Bank observable sky')