16 fits

.github/workflows/build.yml

@@ -5,7 +5,9 @@ on:
   push:
     branches:
     - main
-
+  pull_request:
+    branches:
+    - main
 env:
   BUILD_NAME: 'github-actions[bot]'
   BUILD_EMAIL: '<41898282+github-actions[bot]@users.noreply.github.com>'

build-requirements.txt

@@ -1,3 +1,4 @@
+fitsio
 numpy
 healpy
 camb

examples/advanced/DefineVariables.py

@@ -0,0 +1,90 @@
+'''
+Define Variables
+=======================
+
+To avoid large example files filled with variable definitions we define all relevant variables
+for the advanced exampels here
+
+'''
+
+
+import numpy as np
+import camb
+import glass.ext.camb
+from cosmology import Cosmology
+from glass.fields import generate_lognormal, lognormal_gls
+from glass.lensing import MultiPlaneConvergence
+from glass.shells import distance_grid, tophat_windows
+from glass.observations import gaussian_nz, smail_nz, equal_dens_zbins
+
+
+def get_common_data():
+    global lmax, nside, ncorr, rng, n_arcmin2
+    lmax = nside = 256
+    ncorr = 3
+    # photometric redshift distribution
+    gal_mean_z = [0.5, 1.0]
+    gal_sigma_z = 0.125
+    # total number of galaxies per arcmin2 in each bin
+    gal_dens = 2.0
+    # galaxy number density in units of galaxies/arcmin2/dz
+    z = np.linspace(0., 2., 100)
+    dndz = gaussian_nz(z, gal_mean_z, gal_sigma_z, norm=gal_dens)
+    # ellipticity distribution standard deviation
+    sigma_e = 0.26
+    # photometric redshift error
+    sigma_z0 = 0.03
+    # galaxy bias < 1 so that we can compute expectation from theory
+    beff = 0.9
+    nbin = len(dndz)
+    # labels for tomographic bins
+    binlabels = np.arange(nbin)
+    # name of the catalog we will create
+    catalogName = "myCatalog.FITS"
+    # random number generator, fix seed to make reproducible
+    rng = np.random.default_rng()
+    # galaxy density (using 1/100 of the expected galaxy number density for Stage-IV)
+    n_arcmin2 = 0.3
+    return lmax, nside, sigma_e, sigma_z0, beff, binlabels, catalogName, rng, n_arcmin2
+
+
+def get_glass_data(*, loadCls: bool = False):
+    # cosmology for the simulation
+    h = 0.7
+    Oc = 0.25
+    Ob = 0.05
+    # set up CAMB parameters for matter angular power spectrum
+    pars = camb.set_params(H0=100*h, omch2=Oc*h**2, ombh2=Ob*h**2,
+                           NonLinear=camb.model.NonLinear_both)
+    # get the cosmology from CAMB
+    cosmo = Cosmology.from_camb(pars)
+    # shells of 200 Mpc in comoving distance spacing
+    zgrid = distance_grid(cosmo, 0., 3., dx=200.)
+    # shells from windows
+    shells = tophat_windows(zgrid)
+
+    ws = glass.shells.tophat_windows(zgrid, weight=glass.ext.camb.camb_tophat_weight)
+    if loadCls is True:
+        print("Loading Cls")
+        cls = np.load('../basic/cls.npy')
+    else:
+        print("Calculating Cls (Consider saving and loading them)")
+        cls = glass.ext.camb.matter_cls(pars, lmax, ws)
+
+    gls = lognormal_gls(cls, nside=nside, lmax=lmax, ncorr=ncorr)
+    # assert len(gls) == len(shells) * (len(shells) + 1) // 2
+    # footprint
+
+    print("simulating fields using NSIDE = ", str(nside))
+    # generator for lognormal fields
+    matter = generate_lognormal(gls, nside, ncorr=ncorr, rng=rng)
+    # iterative computation of convergence
+    convergence = MultiPlaneConvergence(cosmo)
+    # true redshift distribution following a Smail distribution
+    z = np.arange(0., 3., 0.01)
+    dndz = smail_nz(z, z_mode=0.9, alpha=2., beta=1.5)
+    dndz *= n_arcmin2
+    # compute tomographic redshift bin edges with equal density
+    nbins = 10
+    zbins = equal_dens_zbins(z, dndz, nbins=nbins)
+    return pars, matter, shells, convergence, ws, z, dndz, zbins

examples/advanced/FITS_IO.py

@@ -0,0 +1,97 @@
+'''
+FITS I/O example
+================
+
+This example creates a galaxy catalogue from a Stage IV Space Satellite Galaxy
+Survey such as *Euclid* and *Roman* combining the :doc:`/basic/plot_density` and
+:doc:`/basic/plot_lensing` examples with galaxy ellipticities and galaxy shears,
+as well as using some auxiliary functions.
+
+The focus in this example is mock catalogue generation using auxiliary functions
+built for simulating Stage IV galaxy surveys.
+
+The catalog is then read and the redshifts plotted
+
+'''
+from itertools import count
+import numpy as np
+from glass.galaxies import galaxy_shear, redshifts_from_nz, gaussian_phz
+from glass.points import positions_from_delta
+from glass.shapes import ellipticity_intnorm
+from glass.lensing import from_convergence
+from glass.user import write_catalog
+from glass.shells import restrict
+from glass.observations import vmap_galactic_ecliptic, tomo_nz_gausserr
+from DefineVariables import get_common_data, get_glass_data
+import healpy as hp
+import matplotlib.pyplot as plt
+import fitsio
+
+
+# Import common variables from create_FITS_helper module
+lmax, nside, sigma_e, sigma_z0, beff, binlabels, catalogName, rng, n_arcmin2 = get_common_data()
+pars, matter, shells, convergence, ws, z, dndz, zbins = get_glass_data(loadCls=False)
+
+# Make a visibility map typical of a space telescope survey, seeing both
+# hemispheres, and low visibility in the galactic and ecliptic bands.
+vis = vmap_galactic_ecliptic(nside)
+# checking the mask:
+hp.mollview(vis, title='Stage IV Space Survey-like Mask', unit='Visibility')
+plt.show()
+
+print("Creating the catalog " + catalogName)
+# open a glass defined write context with a HDU extension name 'CATALOG'
+with write_catalog(catalogName, ext="CATALOG") as out:
+    # go through the shells and simulate
+    for i, delta, shell in zip(count(), matter, shells):
+        z_i, dndz_i = restrict(z, dndz, ws[i])
+        ngal = np.trapz(dndz_i, z_i)
+        print("shell ", i)
+        print("galaxy density: ", ngal)
+        convergence.add_window(delta, shell)
+        kappa = convergence.kappa
+        gamma, = from_convergence(kappa, lmax, shear=True)
+
+        for gal_lon, gal_lat, gal_count in positions_from_delta(ngal, delta, beff, vis):
+            print("galaxies sampled: ", gal_count)
+            gal_eps = ellipticity_intnorm(gal_count, sigma_e, rng=rng)
+            gal_she = galaxy_shear(gal_lon, gal_lat, gal_eps,
+                                   kappa, gamma.real, gamma.imag,
+                                   reduced_shear=False)
+            # some variance in shear weights
+            gal_w = 10**rng.uniform(-2, 2, size=gal_she.shape)
+            gal_z = redshifts_from_nz(gal_count, z_i, dndz_i)
+            # generator photometric redshifts using a Gaussian model
+            gal_phz = gaussian_phz(gal_z, sigma_z0)
+            # attach tomographic bin IDs to galaxies, based on photometric redshifts
+            gal_zbin = np.digitize(gal_phz, np.unique(zbins)) - 1
+
+            out.write(
+                RA=gal_lon,
+                DEC=gal_lat,
+                E1=gal_she.real,
+                E2=gal_she.imag,
+                W=gal_w,
+                Z_TRUE=gal_z,
+                PHZ=gal_phz,
+                ZBIN=gal_zbin,
+            )
+
+# split dndz using the same Gaussian error model assumed in the sampling
+tomo_nz = tomo_nz_gausserr(z, dndz, sigma_z0, zbins)
+nbins = 10
+print("Reading the catalog " + catalogName)
+with fitsio.FITS(catalogName, "r") as fits:
+    data = fits[1].read()
+    plt.figure()
+    plt.title('redshifts in catalogue')
+    plt.ylabel('dN/dz - normalised')
+    plt.xlabel('z')
+    for i in range(nbins):
+        in_bin = (data['ZBIN'] == i)
+        plt.hist(data['Z_TRUE'][in_bin], histtype='stepfilled', edgecolor='none', alpha=0.5, bins=50, density=1, label=f'cat. bin {i}')
+    for i in range(nbins):
+        plt.plot(z, (tomo_nz[i]/n_arcmin2)*nbins, alpha=0.5, label=f'inp. bin {i}')
+    plt.plot(z, dndz/n_arcmin2*nbins, ls='--', c='k')
+    plt.legend(ncol=2)
+    plt.show()