From fb192146923f4a63a3ef93a3fa7a28be378ed33b Mon Sep 17 00:00:00 2001 From: aguinot Date: Tue, 13 Aug 2024 23:43:42 -0400 Subject: [PATCH 1/4] Implement photometry error --- euclidlike_imsim/flux_err.py | 193 +++++++++++++++++++++ euclidlike_imsim/skycat.py | 321 +++++++++++++++++++++++++++++++++++ 2 files changed, 514 insertions(+) create mode 100644 euclidlike_imsim/flux_err.py diff --git a/euclidlike_imsim/flux_err.py b/euclidlike_imsim/flux_err.py new file mode 100644 index 0000000..a583f66 --- /dev/null +++ b/euclidlike_imsim/flux_err.py @@ -0,0 +1,193 @@ +import numpy as np +import numba as nb + +from astropy.time import Time + +import galsim +from galsim.errors import GalSimConfigError + +import ngmix + +import euclidlike + +MIN_IMG_SIZE = 51 +MAX_IMG_SIZE = 501 + +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 + + +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 b940aab..c926b6b 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 @@ -123,6 +126,37 @@ def objects(self): # print('skycat obj 2',process.memory_info().rss) return self._objects + def _build_dtype_dict(self): + self._dtype_dict = {} + obj_types = [] + for coll in self._objects.get_collections(): + objects_type = coll._object_type_unique + if objects_type in obj_types: + continue + col_names = list(coll.native_columns) + for col_name in col_names: + try: + # Some columns cannot be read in snana + np_type = coll.get_native_attribute(col_name).dtype.type() + except Exception as e: + self.logger.warning( + f"The column {col_name} could not be read from skyCatalog." + ) + continue + if np_type is None: + py_type = str + else: + 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. @@ -160,6 +194,177 @@ def getWorldPos(self, index): ra, dec = skycat_obj.ra, skycat_obj.dec return galsim.CelestialCoord(ra * galsim.degrees, dec * galsim.degrees) + def getFlux(self, index, filter=None, mjd=None, exptime=None): + + if filter is None: + filter = self.bandpass.name + if mjd is None: + mjd = self.mjd + if exptime is None: + exptime = self.exptime + + skycat_obj = self.objects[index] + # We cache the SEDs for potential later use + self._seds = skycat_obj.get_observer_sed_components() + for i, sed in enumerate(self._seds.values()): + if i == 0: + sed_sum = sed + else: + sed_sum += sed + raw_flux = skycat_obj.get_euclid_flux( + filter, + sed_sum, + mjd=mjd, + cache=False + ) + if hasattr(skycat_obj, "get_wl_params"): + _, _, mu = skycat_obj.get_wl_params() + else: + mu = 1. + flux = raw_flux * mu * exptime * euclidlike.collecting_area + + 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): + + skycat_obj = self.objects[index] + + if field not in self._dtype_dict: + # We cannot raise an error because one could have a field for snana + # in the config and we don't want to crash because there are no SN + # in this particular image. We then default to False which might not + # be the right type for the required column but we have no way of knowing + # the correct type if the column do not exist. + self.logger.warning(f"The field {field} was not found in skyCatalog.") + return None + elif field not in skycat_obj.native_columns: + if self._dtype_dict[field] is int: + # There are no "special value" for integer so we default to + # hopefully something completely off + return -9999 + elif self._dtype_dict[field] is float: + return np.nan + elif self._dtype_dict[field] is str: + return None + else: + return skycat_obj.get_native_attribute(field) + def getObj(self, index, gsparams=None, rng=None, exptime=30): """ Return the galsim object for the skyCatalog object @@ -345,8 +550,124 @@ def SkyCatWorldPos(config, base, value_type): return pos, safe +def SkyCatValue(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} + params, safe = galsim.config.GetAllParams(config, base, req=req, opt=opt) + field = params["field"] + index = params["index"] + obs_kind = params.get("obs_kind", None) + + if field == "flux": + if obs_kind is None: + val = skycat.getFlux(index) + 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) + val = skycat.getFlux(index, filter=filter, exptime=exptime, mjd=mjd) + else: + val = skycat.getValue(index, field) + + 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") RegisterValueType( "SkyCatWorldPos", SkyCatWorldPos, [galsim.CelestialCoord], input_type="sky_catalog" ) + +# Here we have to provide None as a type otherwise Galsim complains but I don't know why.. +RegisterValueType( + "SkyCatValue", SkyCatValue, [float, int, str, None] # , input_type="sky_catalog" +) +RegisterValueType( + "SkyCatPhotErr", SkyCatPhotErr, [float] # , input_type="sky_catalog" +) From dc520a9e2bf31b3e47a106d95c1b9717b8c22287 Mon Sep 17 00:00:00 2001 From: aguinot Date: Tue, 13 Aug 2024 23:47:39 -0400 Subject: [PATCH 2/4] Update install --- pyproject.toml | 2 ++ 1 file changed, 2 insertions(+) diff --git a/pyproject.toml b/pyproject.toml index a4df064..c8b7834 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", + "numba", ] [project.urls] From d024b550b39abdf90d77721c1bdc30e6af3067f0 Mon Sep 17 00:00:00 2001 From: aguinot Date: Tue, 13 Aug 2024 23:52:09 -0400 Subject: [PATCH 3/4] Fix install --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index c8b7834..1cbc0e2 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -19,7 +19,7 @@ dependencies = [ "numpy>=1.17", # "galsim @ git+https://github.com/GalSim-developers/GalSim#1294_from_images", "astropy>=2.0", - "ngmix", + "ngmix @ git+https://github.com/esheldon/ngmix.git", "numba", ] From 632251e2bebfec29fb143e1cb3e44d8cc8ae86a8 Mon Sep 17 00:00:00 2001 From: aguinot Date: Wed, 21 Aug 2024 12:06:37 -0400 Subject: [PATCH 4/4] Merge main + move nisp info to eulcidlike Add the NISP information to euclidlike `instrument_params.py` file. --- euclidlike/instrument_params.py | 6 +++++ euclidlike_imsim/flux_err.py | 11 ++++----- euclidlike_imsim/skycat.py | 40 +++++++++++++++++++++++++++++++++ 3 files changed, 52 insertions(+), 5 deletions(-) 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 index a583f66..e72172e 100644 --- a/euclidlike_imsim/flux_err.py +++ b/euclidlike_imsim/flux_err.py @@ -5,6 +5,12 @@ 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 @@ -13,11 +19,6 @@ MIN_IMG_SIZE = 51 MAX_IMG_SIZE = 501 -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 - def get_good_img_size(gmix, scale): diff --git a/euclidlike_imsim/skycat.py b/euclidlike_imsim/skycat.py index 0667dc2..2002e6c 100644 --- a/euclidlike_imsim/skycat.py +++ b/euclidlike_imsim/skycat.py @@ -188,6 +188,29 @@ def getWorldPos(self, index): return galsim.CelestialCoord(ra * galsim.degrees, dec * galsim.degrees) def getFlux(self, index, filter=None, mjd=None, exptime=None): + """ + Return the flux associated to an object. + + Parameters + ---------- + index : int + Index of the object in the self.objects catalog. + filter : str, optional + Name of the filter for which the flux is computed. If None, use the + filter provided during initialization. [Default: None] + mjd : float, optional + Date of the observation in MJD format. If None, use the + mjd provided during initialization. [Default: None] + exptime : int or float, optional + Exposure time of the observation. If None, use the + exptime provided during initialization. [Default: None] + + Returns + ------- + flux + Computer flux at the given date for the requested exposure time and + filter. + """ if filter is None: filter = self.bandpass.name @@ -335,6 +358,21 @@ def getFluxErr( return flux_err def getValue(self, index, field): + """ + Return a skyCatalog value for the an object. + + Parameters + ---------- + index : int + Index of the object in the self.objects catalog. + field : str + Name of the field for which you want the value. + + Returns + ------- + int or float or str or None + The value associated to the field or None if the field do not exist. + """ skycat_obj = self.objects[index] @@ -531,6 +569,8 @@ def SkyCatWorldPos(config, base, value_type): def SkyCatValue(config, base, value_type): + """Return a value from the object part of the skyCatalog + """ skycat = galsim.config.GetInputObj("sky_catalog", config, base, "SkyCatValue")