How to use astropy - 10 common examples

'-o', '--output', metavar='OUTPUT.fits[.gz]', required=True,
    help='name of output FITS file [required]')
opts = parser.parse_args()


# Late imports.
import healpy as hp
import numpy as np
from astropy.table import Table
from lalinference.io import fits
from lalinference.bayestar import distance
from lalinference.bayestar import moc
from lalinference.bayestar.sky_map import derasterize
from lalinference.healpix_tree import adaptive_healpix_histogram

samples = Table.read(opts.input, format='ascii')
theta = 0.5*np.pi - samples['dec']
phi = samples['ra']
if 'distance' in samples.colnames:
    samples.rename_column('distance', 'dist')

p = adaptive_healpix_histogram(
    theta, phi, opts.samples_per_bin,
    nside=opts.nside, max_nside=opts.max_nside, nest=True)


def diststats(samples, max_nside, nside, ipix):
    step = (max_nside // nside) ** 2
    i0 = np.searchsorted(samples['ipix'], step * ipix)
    i1 = np.searchsorted(samples['ipix'], step * (ipix + 1))
    if i0 == i1:
        return np.inf, 0.0

if isinstance(par.value, numbers.Number):
                ov = par.value
                if isinstance(par, mp.MJDParameter):
                    continue
                else:
                    par.value = ov * 0.8
        self.res = Residuals(
            self.toasB1855, self.modelB1855, use_weighted_mean=False
        ).time_resids.to(u.s)
        f = open(self.out_parfile, "w")
        f.write(self.modelB1855.as_parfile())
        f.close()
        read_model = mb.get_model(self.out_parfile)
        read_res = Residuals(
            self.toasB1855, read_model, use_weighted_mean=False
        ).time_resids.to(u.s)
        assert np.all(
            np.abs(read_res.value - self.res.value) < 1e-15
        ), "Output parfile did not produce same residuals."
        for pp in self.modelB1855.params:
            par_ori = getattr(self.modelB1855, pp)
            par_read = getattr(read_model, pp)
            if par_ori.uncertainty_value is not None:
                unc_diff = par_ori.uncertainty_value - par_read.uncertainty_value
                assert np.abs(unc_diff) < 1e-15, (
                    pp
                    + "uncertainty does not keep the precision. at"
                    + str(np.abs(unc_diff))
                )

def test_vela(self):
        print("Test RMS of a VELA ephemeris with WAVE parameters.")
        rs = pint.residuals.Residuals(self.t, self.m).time_resids
        rms = rs.to(u.us).std()
        emsg = "RMS of " + str(rms.value) + " is too big."
        assert rms < 350.0 * u.us, emsg

def test_fill_cube():
    filename = gammapy_extra.filename('test_datasets/background/bg_cube_model_test1.fits')
    array = Cube.read(filename, format='table', scheme='bg_cube')
    array.data = Quantity(np.zeros_like(array.data.value), 'u')
    print(type(array.data))

    dir = str(gammapy_extra.dir) + '/datasets/hess-crab4-hd-hap-prod2'
    data_store = DataStore.from_dir(dir)
    ev_list = data_store.load_all('events')

    array.fill_events(ev_list)

    array.write('test_background.fits', format='image', clobber=True)

array = np.array([[1.0, 1.0, 1.0], [1.0, 2.0, 1.0], [1.0, 1.0, 1.0]])

fits.writeto(data=array, filename=path + "3x3_ones_central_two.fits")
stop

array1 = np.ones((3, 3))
array2 = 2.0 * np.ones((3, 3))
array3 = 3.0 * np.ones((3, 3))
array4 = 4.0 * np.ones((3, 3))
array5 = 5.0 * np.ones((3, 3))
array6 = 6.0 * np.ones((3, 3))
array7 = 7.0 * np.ones((3, 3))
array8 = 8.0 * np.ones((3, 3))

fits.writeto(data=array2, filename=path + "3x3_twos.fits")
fits.writeto(data=array3, filename=path + "3x3_threes.fits")
fits.writeto(data=array4, filename=path + "3x3_fours.fits")
fits.writeto(data=array5, filename=path + "3x3_fives.fits")
fits.writeto(data=array6, filename=path + "3x3_sixes.fits")
fits.writeto(data=array7, filename=path + "3x3_sevens.fits")
fits.writeto(data=array8, filename=path + "3x3_eights.fits")

new_hdul = fits.HDUList()
new_hdul.append(fits.ImageHDU(array1))
new_hdul.append(fits.ImageHDU(array2))
new_hdul.append(fits.ImageHDU(array3))
new_hdul.append(fits.ImageHDU(array4))
new_hdul.append(fits.ImageHDU(array5))
new_hdul.append(fits.ImageHDU(array6))
new_hdul.append(fits.ImageHDU(array7))
new_hdul.append(fits.ImageHDU(array8))

def astropySkyToPixels(self, astropyWcs, skyPosList):
        """Use an astropy wcs to convert pixels to sky

        @param[in] astropyWcs  a celestial astropy.wcs.WCS with 2 axes in RA, Dec order
        @param[in] skyPosList ICRS sky coordinates as a list of lsst.geom.SpherePoint
        @returns a list of lsst.geom.Point2D, 0-based pixel positions

        Converts the input from ICRS to the coordinate system of the wcs
        """
        skyCoordList = [astropy.coordinates.SkyCoord(c[0].asDegrees(),
                                                     c[1].asDegrees(),
                                                     frame="icrs",
                                                     unit="deg") for c in skyPosList]
        xyArr = [astropy.wcs.utils.skycoord_to_pixel(coords=sc,
                                                     wcs=astropyWcs,
                                                     origin=0,
                                                     mode="all") for sc in skyCoordList]
        # float is needed to avoid truncation to int
        return [lsst.geom.Point2D(float(x), float(y)) for x, y in xyArr]

def tmpFits(*hdus):
    # Given a list of numpy arrays, create a temporary FITS file that
    # contains them as consecutive HDUs. Yield it, then remove it.
    hdus = [astropy.io.fits.PrimaryHDU(hdus[0])] + \
        [astropy.io.fits.ImageHDU(hdu) for hdu in hdus[1:]]
    hdulist = astropy.io.fits.HDUList(hdus)
    tempdir = tempfile.mkdtemp()
    try:
        filename = os.path.join(tempdir, 'test.fits')
        hdulist.writeto(filename)
        yield filename
    finally:
        shutil.rmtree(tempdir)

self.assertAlmostEqual(angle.EL,el_iod4.deg,6, msg="IOD Format 4 failed in EL")
        self.assertEqual(angle.Epoch,2000)

        # 		 5: AZ/EL  = DDDMMmm+DDMMmm MX   (MX in minutes of arc)
        #              IOD 5 1835170-081148 (1855 epoch)
        az_iod5  = Angle('183d51.70m')  # (DDDMMmm)
        el_iod5 = Angle('-08d11.48m') # (DDMMmm)

        angle = iod.Angle(5,5,'1835170-081148',0,"","IOD")
        self.assertAlmostEqual(angle.AZ,az_iod5.deg,6, msg="IOD Format 5 failed in AZ")
        self.assertAlmostEqual(angle.EL,el_iod5.deg,6, msg="IOD Format 5 failed in EL")
        self.assertEqual(angle.Epoch,2000)

        # 		 6: AZ/EL  = DDDdddd+DDdddd MX   (MX in degrees of arc)
        #              IOD 6 0214030+732900
        az_iod6  = Angle('021.4030d')  # (DDDdddd)
        el_iod6 = Angle('73.2900d')   # (DDdddd)

        angle = iod.Angle(6,5,'0214030+732900',0,"","IOD")
        self.assertAlmostEqual(angle.AZ,az_iod6.deg,6, msg="IOD Format 6 failed in AZ")
        self.assertAlmostEqual(angle.EL,el_iod6.deg,6, msg="IOD Format 6 failed in EL")
        self.assertEqual(angle.Epoch,2000)

        # 		 7: RA/DEC = HHMMSSs+DDdddd MX   (MX in degrees of arc)
        #              IOD 7 2047449+293762
        ra_iod7  = Angle('20h47m44.9s')  # (HHMMSSs)
        dec_iod7 = Angle('29.3762d')     # (DDdddd)

        angle = iod.Angle(7,5,'2047449+293762',0,"","IOD")
        self.assertAlmostEqual(angle.RA,ra_iod7.deg,6, msg="IOD Format 7 failed in RA")
        self.assertAlmostEqual(angle.DEC,dec_iod7.deg,6, msg="IOD Format 7 failed in DEC")
        self.assertEqual(angle.Epoch,2000)

"orbit": [
        26553.4 * u.km,
        0.741 * u.one,
        63.4 * u.deg,
        0.0 * u.deg,
        -10.12921 * u.deg,
        0.0 * u.rad,
    ],
    "period": 28 * u.day,
}

moon_leo = {
    "body": Moon,
    "tof": 60 * u.day,
    "raan": -2.18 * 1e-4 * u.deg,
    "argp": 15.0 * 1e-3 * u.deg,
    "inc": 6.0 * 1e-4 * u.deg,
    "orbit": [
        6678.126 * u.km,
        0.01 * u.one,
        28.5 * u.deg,
        0.0 * u.deg,
        0.0 * u.deg,
        0.0 * u.rad,
    ],
    "period": 28 * u.day,
}

moon_geo = {
    "body": Moon,
    "tof": 60 * u.day,
    "raan": 6.0 * u.deg,

did_tempo1 = False
    try:
        import tempo_utils

        log.info("Running TEMPO1...")
        t1_result = np.genfromtxt(
            t1_parfile + ".tempo_test", names=True, comments="#", dtype=np.longdouble
        )
        t1_resids = t1_result["residuals_phase"] / float(m.F0.value) * 1e6 * u.us
        did_tempo1 = True
        diff_t1 = (resids_us - t1_resids).to(u.ns)
        diff_t1 -= diff_t1.mean()
        log.info(
            "Max resid diff between PINT and T1: %.2f ns" % np.fabs(diff_t1).max().value
        )
        log.info("Std resid diff between PINT and T1: %.2f ns" % diff_t1.std().value)
        diff_t2_t1 = (t2_resids - t1_resids).to(u.ns)
        diff_t2_t1 -= diff_t2_t1.mean()
        log.info(
            "Max resid diff between T1 and T2: %.2f ns"
            % np.fabs(diff_t2_t1).max().value
        )
        log.info("Std resid diff between T1 and T2: %.2f ns" % diff_t2_t1.std().value)
    except:
        pass

    if did_tempo1 and not planets:
        assert np.fabs(diff_t1).max() < 32.0 * u.ns

    def do_plot():
        plt.clf()
        plt.subplot(211)