diff --git a/euclidlike/instrument_params.py b/euclidlike/instrument_params.py index b24abfe..7c4a46d 100644 --- a/euclidlike/instrument_params.py +++ b/euclidlike/instrument_params.py @@ -115,6 +115,12 @@ vis_bands = ['VIS'] nisp_bands = ['NISP_Y', 'NISP_J', 'NISP_H'] +# Add some information related to NISP +nisp_gain = 2 # https://arxiv.org/pdf/2405.13496 Sect 4.3.7 +nisp_pixel_scale = 0.3 +nisp_dark_current = 0.02 # e.s-1.pix-1 https://arxiv.org/pdf/2405.13493 Sect 4.1.2 +nisp_read_noise = 6.2 # e.pix-1 https://arxiv.org/pdf/2405.13496 Sect 4.3.5 + # Items to potentially do later; part of the galsim.roman setup that currently has no correspondence # here. # dark_current diff --git a/euclidlike_imsim/flux_err.py b/euclidlike_imsim/flux_err.py new file mode 100644 index 0000000..e72172e --- /dev/null +++ b/euclidlike_imsim/flux_err.py @@ -0,0 +1,194 @@ +import numpy as np +import numba as nb + +from astropy.time import Time + +import galsim +from galsim.errors import GalSimConfigError +from euclidlike.instrument_params import ( + nisp_gain, + nisp_pixel_scale, + nisp_dark_current, + nisp_read_noise, +) + +import ngmix + +import euclidlike + +MIN_IMG_SIZE = 51 +MAX_IMG_SIZE = 501 + + +def get_good_img_size(gmix, scale): + + _, _, sigma = gmix.get_e1e2sigma() + size_pix = sigma / scale + img_size = min(max(size_pix * 5, MIN_IMG_SIZE), MAX_IMG_SIZE) + img_size = int(np.ceil(img_size)) + return img_size + + +@nb.njit +def get_flux(gmix, pixels): + aper_flux = 0 + n_pixels = pixels.shape[0] + for i in range(n_pixels): + aper_flux += ngmix.gmix.gmix_nb.gmix_eval_pixel_fast(gmix, pixels[i]) + return aper_flux + + +@nb.njit +def circular_aper(N, r): + """ + Draws a circle of radius r in a square image of size N x N. + + Parameters: + N (int): The size of the square image (N x N). + r (int): The radius of the circle. + + Returns: + numpy.ndarray: A 2D numpy array representing the image with the circle. + """ + # Create an empty N x N array (filled with zeros) + image = np.zeros((N, N), dtype=np.int16) + + # Define the center of the image + center = N // 2 + + # Iterate over each pixel in the image + for y in range(N): + for x in range(N): + # Calculate the distance from the center + distance = np.sqrt((x - center)**2 + (y - center)**2) + + # If the distance is less than or equal to the radius, set the pixel to 1 + if distance <= r: + image[y, x] = 1 + + return image + + +def get_variance(bandpass, world_pos, mjd, exptime): + + filter = bandpass.name + + sky_lvl = euclidlike.getSkyLevel( + bandpass, + world_pos=world_pos, + date=Time(mjd, format="mjd").datetime, + exptime=exptime, + ) + + if filter == "VIS": + return ( + sky_lvl * euclidlike.pixel_scale**2 + + euclidlike.read_noise**2 + ) / euclidlike.gain**2 + elif "NISP" in filter: + return ( + sky_lvl * nisp_pixel_scale**2 + + nisp_read_noise**2 + + (nisp_dark_current * exptime)**2 + ) / nisp_gain**2 + + +def get_flux_err( + ra, + dec, + wcs, + bandpass, + mjd, + exptime, + aper_type, + aper_size, + Nexp, + gal_pars=None, + star_flux=None, + model=None, +): + if ( + (gal_pars is None and star_flux is None) + | (gal_pars is not None and star_flux is not None) + ): + raise ValueError( + "One of [gal_pars, star_flux] must be provided, not both." + ) + if model is None: + raise ValueError("Model must be provided in ['bd', 'gausse']") + + filter = bandpass.name + wave = bandpass.effective_wavelength + + lam_over_diam = np.rad2deg( + wave * 1e-9 / euclidlike.diameter + ) * 3600 + psf_pars = [ + 0, 0, 0, 0, 2 * (0.5 * lam_over_diam)**2, 1 + ] + + if gal_pars is not None: + gal_gmix = ngmix.gmix.make_gmix_model(gal_pars, model) + psf_gmix = ngmix.gmix.make_gmix_model(psf_pars, "gauss") + obj_gmix = gal_gmix.convolve(psf_gmix) + else: + point_source_pars = psf_pars + point_source_pars[-1] = star_flux + obj_gmix = ngmix.gmix.make_gmix_model(point_source_pars, "gauss") + + world_pos = galsim.CelestialCoord( + ra=ra * galsim.degrees, dec=dec * galsim.degrees + ) + if filter == "VIS": + jacobian = wcs.jacobian(world_pos=world_pos) + else: + jacobian = wcs.jacobian() + pixel_scale = np.sqrt(jacobian.pixelArea()) + + img_size = get_good_img_size(obj_gmix, pixel_scale) + + jacob = ngmix.Jacobian( + row=(img_size - 1) / 2, + col=(img_size - 1) / 2, + wcs=jacobian, + ) + + gmix_img = obj_gmix.make_image( + (img_size, img_size), + jacobian=jacob, + fast_exp=True, + ) + + noise_var = get_variance(bandpass, world_pos, mjd, exptime) + if aper_type == "circular": + if aper_size < 0: + aperture_mask = np.ones((img_size, img_size)) + else: + aperture_mask = circular_aper(img_size, aper_size / pixel_scale) + else: + raise NotImplementedError + weight_img = ( + aperture_mask + * np.ones_like(gmix_img) / ((noise_var / Nexp + gmix_img)) + ) + pixels = ngmix.pixels.make_pixels( + gmix_img, + weight_img, + jacob, + ) + + flux = get_flux( + obj_gmix.get_data(), + pixels, + ) + + snr = np.sqrt( + ngmix.gmix.gmix_nb.get_model_s2n_sum( + obj_gmix.get_data(), + pixels, + ) + ) + + flux_err = flux / snr + + return flux, flux_err, snr diff --git a/euclidlike_imsim/skycat.py b/euclidlike_imsim/skycat.py index 9474036..2002e6c 100644 --- a/euclidlike_imsim/skycat.py +++ b/euclidlike_imsim/skycat.py @@ -11,6 +11,9 @@ RegisterValueType, RegisterObjectType, ) +from galsim.errors import GalSimConfigValueError + +from .flux_err import get_flux_err import euclidlike @@ -139,6 +142,14 @@ def _build_dtype_dict(self): py_type = type(np_type.astype(object)) self._dtype_dict[col_name] = py_type + def _get_bandpass(self, filter): + + if hasattr(self, "_bandpass"): + return self._bandpass[filter] + else: + self._bandpass = euclidlike.getBandpasses() + return self._bandpass[filter] + def get_ccd_center(self): """ Return the CCD center. @@ -230,6 +241,122 @@ def getFlux(self, index, filter=None, mjd=None, exptime=None): return flux + def getFluxErr( + self, + index, + filter=None, + mjd=None, + exptime=None, + aper_type="circular", + aper_size=1, + Nexp=1, + ): + + if filter is None: + filter = self.bandpass.name + bandpass = self.bandpass + else: + bandpass = self._get_bandpass(filter) + if mjd is None: + mjd = self.mjd + if exptime is None: + exptime = self.exptime + + skycat_obj = self.objects[index] + + # seds = skycat_obj.get_observer_sed_components() + seds = self._seds + + if filter == "VIS": + gain = euclidlike.gain + wcs = self.wcs + else: + gain = 2 + exptime = 87.2 + wcs = galsim.PixelScale(scale=0.3) + + if skycat_obj.object_type == "diffsky_galaxy": + _, _, mu = skycat_obj.get_wl_params() + disk_hlr = skycat_obj.get_native_attribute("diskHalfLightRadiusArcsec") + bulge_hlr = skycat_obj.get_native_attribute("spheroidHalfLightRadiusArcsec") + disk_T = 2 * (disk_hlr / 1.17741)**2 + bulge_T = 2 * (bulge_hlr / 1.17741)**2 + flux_disk = ( + skycat_obj.get_euclid_flux( + filter, + sed=seds["disk"] + seds["knots"], + mjd=mjd, + cache=False + ) + * mu + * exptime + * euclidlike.collecting_area + / gain + ) + flux_bulge = ( + skycat_obj.get_euclid_flux( + filter, + seds["bulge"], + mjd=mjd, + cache=False + ) + * mu + * exptime + * euclidlike.collecting_area + / gain + ) + flux_tot = flux_disk + flux_bulge + g1 = skycat_obj.get_native_attribute("diskEllipticity1") + g2 = skycat_obj.get_native_attribute("diskEllipticity2") + T = disk_T + gal_pars = [ + 0, # center row + 0, # center col + g1, # shear g1 + g2, # shear g2 + T, # total size + np.log10(bulge_T / disk_T), # log10(T_dev/T_exp) + flux_bulge / flux_tot, # fluxfracdev F_dev/F_tot + flux_tot, # total flux + ] + star_flux = None + model = "bd" + else: + flux_tot = ( + skycat_obj.get_euclid_flux( + filter, + sed=seds["this_object"], + mjd=mjd, + cache=False, + ) + * exptime + * euclidlike.collecting_area + / gain + ) + gal_pars = None + star_flux = flux_tot + model = "gauss" + + ra = skycat_obj.ra + dec = skycat_obj.dec + + flux_err = get_flux_err( + ra=ra, + dec=dec, + wcs=wcs, + bandpass=bandpass, + mjd=mjd, + exptime=exptime, + aper_type=aper_type, + aper_size=aper_size, + Nexp=Nexp, + gal_pars=gal_pars, + star_flux=star_flux, + model=model, + ) + + return flux_err + def getValue(self, index, field): """ Return a skyCatalog value for the an object. @@ -475,6 +602,81 @@ def SkyCatValue(config, base, value_type): return val, safe +def SkyCatPhotErr(config, base, value_type): + + skycat = galsim.config.GetInputObj("sky_catalog", config, base, "SkyCatValue") + + # Setup the indexing sequence if it hasn't been specified. The + # normal thing with a catalog is to just use each object in order, + # so we don't require the user to specify that by hand. We can do + # it for them. + galsim.config.SetDefaultIndex(config, skycat.getNObjects()) + + req = {"field": str, "index": int} + opt = {"obs_kind": str, "aper_type": str, "aper_size": float, "Nexp": int} + params, safe = galsim.config.GetAllParams(config, base, req=req, opt=opt) + field = params["field"] + index = params["index"] + obs_kind = params.get("obs_kind", None) + aper_type = params.get("aper_type", "circular") + aper_size = params.get("aper_size", 1.) + Nexp = params.get("Nexp", 1) + + if field not in ["flux", "flux_err", "snr"]: + raise GalSimConfigValueError( + "Invalid SkyCatPhotErr field.", field, ["flux", "flux_err", "snr"] + ) + + # Check for cached values + if "phot_err" in base: + if base["obj_num"] == base["phot_err"]["obj_num"]: + phot_err = base["phot_err"] + # We check the config is the same + if ( + (phot_err["obs_kind"] == obs_kind) + & (phot_err["aper_type"] == aper_type) + & (phot_err["aper_size"] == aper_size) + & (phot_err["Nexp"] == Nexp) + ): + return phot_err[field] + + if obs_kind is None: + filter = None + exptime = None + mjd = None + else: + pointing = galsim.config.GetInputObj( + "obseq_data", config, base, "OpSeqDataLoader" + ) + filter = pointing.get("filter", obs_kind=obs_kind) + exptime = pointing.get("exptime", obs_kind=obs_kind) + mjd = pointing.get("mjd", obs_kind=obs_kind) + + flux, flux_err, snr = skycat.getFluxErr( + index, + filter=filter, + exptime=exptime, + mjd=mjd, + aper_type=aper_type, + aper_size=aper_size, + Nexp=Nexp, + ) + + # Cache value + phot_err = { + "obj_num": base["obj_num"], + "obs_kind": obs_kind, + "aper_type": aper_type, + "aper_size": aper_size, + "Nexp": Nexp, + "flux": flux, + "flux_err": flux_err, + "snr": snr, + } + base["phot_err"] = phot_err + + return phot_err[field], safe + RegisterInputType("sky_catalog", SkyCatalogLoader(SkyCatalogInterface, has_nobj=True)) RegisterObjectType("SkyCatObj", SkyCatObj, input_type="sky_catalog") @@ -486,3 +688,6 @@ def SkyCatValue(config, base, value_type): RegisterValueType( "SkyCatValue", SkyCatValue, [float, int, str, None] # , input_type="sky_catalog" ) +RegisterValueType( + "SkyCatPhotErr", SkyCatPhotErr, [float] # , input_type="sky_catalog" +) diff --git a/pyproject.toml b/pyproject.toml index a4df064..1cbc0e2 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -19,6 +19,8 @@ dependencies = [ "numpy>=1.17", # "galsim @ git+https://github.com/GalSim-developers/GalSim#1294_from_images", "astropy>=2.0", + "ngmix @ git+https://github.com/esheldon/ngmix.git", + "numba", ] [project.urls]