From efd8a67b7a1cdb43b380e9c13bc4f722050583c2 Mon Sep 17 00:00:00 2001 From: Eli Rykoff Date: Thu, 19 Oct 2023 19:13:36 -0700 Subject: [PATCH 01/11] Add DES to transformed catalogs. --- python/lsst/the/monster/isolate_and_transform.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/python/lsst/the/monster/isolate_and_transform.py b/python/lsst/the/monster/isolate_and_transform.py index 0292bf7..c0739db 100644 --- a/python/lsst/the/monster/isolate_and_transform.py +++ b/python/lsst/the/monster/isolate_and_transform.py @@ -4,7 +4,7 @@ from lsst.pipe.tasks.isolatedStarAssociation import IsolatedStarAssociationTask from .splinecolorterms import ColortermSpline -from .refcats import GaiaXPInfo, GaiaDR3Info, SkyMapperInfo, PS1Info, VSTInfo +from .refcats import GaiaXPInfo, GaiaDR3Info, SkyMapperInfo, PS1Info, VSTInfo, DESInfo from .utils import read_stars, makeRefSchema, makeRefCat __all__ = ["MatchAndTransform"] @@ -42,7 +42,7 @@ class MatchAndTransform: def __init__(self, gaia_reference_class=GaiaDR3Info, catalog_info_class_list=[GaiaXPInfo, SkyMapperInfo, - PS1Info, VSTInfo], + PS1Info, VSTInfo, DESInfo], write_path_inp=None, testing_mode=False, ): From f03eb5fbad6973e86f37c9a2e4fa158a9f38565e Mon Sep 17 00:00:00 2001 From: Eli Rykoff Date: Thu, 19 Oct 2023 19:14:09 -0700 Subject: [PATCH 02/11] Return empty list if no stars are found. --- python/lsst/the/monster/utils.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/python/lsst/the/monster/utils.py b/python/lsst/the/monster/utils.py index 0377ea8..89d9d70 100644 --- a/python/lsst/the/monster/utils.py +++ b/python/lsst/the/monster/utils.py @@ -43,6 +43,9 @@ def read_stars(path, indices, allow_missing=False): else: stars.extend(temp) + if stars is None: + return [] + stars = stars.copy(deep=True).asAstropy() stars["coord_ra"].convert_unit_to(units.degree) From 2ff62a70bcd61bf059fe590e2e194c118b208217 Mon Sep 17 00:00:00 2001 From: Eli Rykoff Date: Thu, 19 Oct 2023 19:15:20 -0700 Subject: [PATCH 03/11] Add PS1 to DES colorterms to test directory. --- .../colorterms/TestPS1_to_DES_band_g.yaml | 50 +++++++++++++++++++ .../colorterms/TestPS1_to_DES_band_i.yaml | 50 +++++++++++++++++++ .../colorterms/TestPS1_to_DES_band_r.yaml | 50 +++++++++++++++++++ .../colorterms/TestPS1_to_DES_band_y.yaml | 50 +++++++++++++++++++ .../colorterms/TestPS1_to_DES_band_z.yaml | 50 +++++++++++++++++++ tests/test_catalog_measurement.py | 1 + 6 files changed, 251 insertions(+) create mode 100644 tests/data/colorterms/TestPS1_to_DES_band_g.yaml create mode 100644 tests/data/colorterms/TestPS1_to_DES_band_i.yaml create mode 100644 tests/data/colorterms/TestPS1_to_DES_band_r.yaml create mode 100644 tests/data/colorterms/TestPS1_to_DES_band_y.yaml create mode 100644 tests/data/colorterms/TestPS1_to_DES_band_z.yaml diff --git a/tests/data/colorterms/TestPS1_to_DES_band_g.yaml b/tests/data/colorterms/TestPS1_to_DES_band_g.yaml new file mode 100644 index 0000000..9ddefff --- /dev/null +++ b/tests/data/colorterms/TestPS1_to_DES_band_g.yaml @@ -0,0 +1,50 @@ +flux_offset: 9.383807258485742 +mag_nodes: +- 13.0 +- 13.7 +- 14.4 +- 15.1 +- 15.8 +- 16.5 +- 17.2 +- 17.9 +- 18.6 +- 19.3 +mag_spline_values: +- 1.021093776662617 +- 1.01613939998474 +- 1.0153732139674558 +- 1.0134087220962626 +- 1.0102851182068038 +- 1.0074847295795457 +- 1.003988315134559 +- 1.0003494218256264 +- 1.001014866965003 +- 0.9953993922059542 +nodes: +- 0.5 +- 0.7555555555555555 +- 1.011111111111111 +- 1.2666666666666666 +- 1.5222222222222221 +- 1.7777777777777777 +- 2.033333333333333 +- 2.2888888888888888 +- 2.5444444444444443 +- 2.8 +source_color_field_1: g_flux +source_color_field_2: i_flux +source_field: g_flux +source_survey: PS1 +spline_values: +- 0.9751839369285394 +- 0.9629421889103474 +- 0.9521908858321729 +- 0.9434778630350008 +- 0.935988590058585 +- 0.9310037854643703 +- 0.9274918780814252 +- 0.92806696150321 +- 0.9281549226667183 +- 0.9314907727679823 +target_survey: DES diff --git a/tests/data/colorterms/TestPS1_to_DES_band_i.yaml b/tests/data/colorterms/TestPS1_to_DES_band_i.yaml new file mode 100644 index 0000000..7509a06 --- /dev/null +++ b/tests/data/colorterms/TestPS1_to_DES_band_i.yaml @@ -0,0 +1,50 @@ +flux_offset: -71.36179718638563 +mag_nodes: +- 13.75 +- 14.166666666666666 +- 14.583333333333334 +- 15.0 +- 15.416666666666666 +- 15.833333333333334 +- 16.25 +- 16.666666666666668 +- 17.083333333333332 +- 17.5 +mag_spline_values: +- 1.0113305045934213 +- 1.0060211085345212 +- 1.0059145859254501 +- 1.0058228250081558 +- 1.0042729735075449 +- 1.0029713243411502 +- 1.001076395498808 +- 0.9989529197851542 +- 0.9977368461225616 +- 0.9971579189496544 +nodes: +- 0.5 +- 0.7555555555555555 +- 1.011111111111111 +- 1.2666666666666666 +- 1.5222222222222221 +- 1.7777777777777777 +- 2.033333333333333 +- 2.2888888888888888 +- 2.5444444444444443 +- 2.8 +source_color_field_1: g_flux +source_color_field_2: i_flux +source_field: i_flux +source_survey: PS1 +spline_values: +- 1.0005976475311547 +- 1.0122185576327776 +- 1.0233149953987712 +- 1.036603558419902 +- 1.0536876004201725 +- 1.076813267264628 +- 1.1094900370659144 +- 1.149170115645006 +- 1.187223154498976 +- 1.2229387496358388 +target_survey: DES diff --git a/tests/data/colorterms/TestPS1_to_DES_band_r.yaml b/tests/data/colorterms/TestPS1_to_DES_band_r.yaml new file mode 100644 index 0000000..fe1d1e4 --- /dev/null +++ b/tests/data/colorterms/TestPS1_to_DES_band_r.yaml @@ -0,0 +1,50 @@ +flux_offset: -2.618738793757136 +mag_nodes: +- 13.25 +- 13.777777777777779 +- 14.305555555555555 +- 14.833333333333334 +- 15.36111111111111 +- 15.88888888888889 +- 16.416666666666668 +- 16.944444444444443 +- 17.47222222222222 +- 18.0 +mag_spline_values: +- 1.0248123394811852 +- 1.0132226522734895 +- 1.0104573485590365 +- 1.0095428985646504 +- 1.0075950739037738 +- 1.0053214808492152 +- 1.0022975482439591 +- 0.999771151064734 +- 0.9983252191675076 +- 0.9968004690652866 +nodes: +- 0.5 +- 0.7555555555555555 +- 1.011111111111111 +- 1.2666666666666666 +- 1.5222222222222221 +- 1.7777777777777777 +- 2.033333333333333 +- 2.2888888888888888 +- 2.5444444444444443 +- 2.8 +source_color_field_1: g_flux +source_color_field_2: i_flux +source_field: r_flux +source_survey: PS1 +spline_values: +- 1.0202994425388818 +- 1.0347397306628556 +- 1.0487873695700998 +- 1.062591470599521 +- 1.0783049215980738 +- 1.096656958217486 +- 1.1209192430189703 +- 1.1537085199184196 +- 1.1915634746104733 +- 1.2280472756809007 +target_survey: DES diff --git a/tests/data/colorterms/TestPS1_to_DES_band_y.yaml b/tests/data/colorterms/TestPS1_to_DES_band_y.yaml new file mode 100644 index 0000000..f5c7e01 --- /dev/null +++ b/tests/data/colorterms/TestPS1_to_DES_band_y.yaml @@ -0,0 +1,50 @@ +flux_offset: -277.47544114634593 +mag_nodes: +- 13.75 +- 14.11111111111111 +- 14.472222222222221 +- 14.833333333333334 +- 15.194444444444445 +- 15.555555555555555 +- 15.916666666666666 +- 16.27777777777778 +- 16.63888888888889 +- 17.0 +mag_spline_values: +- 0.9848489057200971 +- 0.9850494710367231 +- 0.9863914427836119 +- 0.988381958233075 +- 0.991931260351768 +- 0.9949739344290779 +- 0.9994894266911729 +- 1.0003893893671911 +- 1.0034863438619643 +- 0.9992251418647077 +nodes: +- 0.0 +- 0.07777777777777778 +- 0.15555555555555556 +- 0.23333333333333334 +- 0.3111111111111111 +- 0.3888888888888889 +- 0.4666666666666667 +- 0.5444444444444445 +- 0.6222222222222222 +- 0.7 +source_color_field_1: i_flux +source_color_field_2: z_flux +source_field: y_flux +source_survey: PS1 +spline_values: +- 0.9599292090263586 +- 0.9702373734684031 +- 0.9788308188533541 +- 0.9821375908922485 +- 0.985943098794281 +- 0.9890395284463673 +- 0.9927616682895601 +- 0.9975948708557999 +- 1.002450609732945 +- 1.0085672495501865 +target_survey: DES diff --git a/tests/data/colorterms/TestPS1_to_DES_band_z.yaml b/tests/data/colorterms/TestPS1_to_DES_band_z.yaml new file mode 100644 index 0000000..b96fb6a --- /dev/null +++ b/tests/data/colorterms/TestPS1_to_DES_band_z.yaml @@ -0,0 +1,50 @@ +flux_offset: -23.665562628998543 +mag_nodes: +- 13.75 +- 14.166666666666666 +- 14.583333333333334 +- 15.0 +- 15.416666666666666 +- 15.833333333333334 +- 16.25 +- 16.666666666666668 +- 17.083333333333332 +- 17.5 +mag_spline_values: +- 1.0049538428509763 +- 1.0055001925777558 +- 1.0054234794424295 +- 1.0036491033255412 +- 1.0023833924181555 +- 1.0010346171147555 +- 0.9999443463804955 +- 0.9987417125909357 +- 0.9977394356884955 +- 0.9976415031801529 +nodes: +- 0.0 +- 0.07777777777777778 +- 0.15555555555555556 +- 0.23333333333333334 +- 0.3111111111111111 +- 0.3888888888888889 +- 0.4666666666666667 +- 0.5444444444444445 +- 0.6222222222222222 +- 0.7 +source_color_field_1: i_flux +source_color_field_2: z_flux +source_field: z_flux +source_survey: PS1 +spline_values: +- 1.0096298727782347 +- 1.0277244564625445 +- 1.053427229781647 +- 1.070013998965467 +- 1.085958750417689 +- 1.0999856044852463 +- 1.1160373581065222 +- 1.1339253457491285 +- 1.1557951134586582 +- 1.1778150058178611 +target_survey: DES diff --git a/tests/test_catalog_measurement.py b/tests/test_catalog_measurement.py index 49a2867..16771d3 100644 --- a/tests/test_catalog_measurement.py +++ b/tests/test_catalog_measurement.py @@ -34,6 +34,7 @@ class DESInfoTester(DESInfo): class PS1InfoTester(PS1Info): PATH = os.path.join(ROOT, "data", "ps1") NAME = "TestPS1" + COLORTERM_PATH = os.path.join(ROOT, "data", "colorterms") class GaiaXPSplineMeasurerTester(GaiaXPSplineMeasurer): From 3f35ac892bfc640e9bea088ddd80b4eadc91cc12 Mon Sep 17 00:00:00 2001 From: Eli Rykoff Date: Thu, 19 Oct 2023 19:15:40 -0700 Subject: [PATCH 04/11] Add code to make and plot offset maps. --- bin.src/make_all_offset_maps.py | 27 ++ python/lsst/the/monster/__init__.py | 1 + python/lsst/the/monster/measure_offsetmaps.py | 321 ++++++++++++++++++ 3 files changed, 349 insertions(+) create mode 100644 bin.src/make_all_offset_maps.py create mode 100644 python/lsst/the/monster/measure_offsetmaps.py diff --git a/bin.src/make_all_offset_maps.py b/bin.src/make_all_offset_maps.py new file mode 100644 index 0000000..6c0a406 --- /dev/null +++ b/bin.src/make_all_offset_maps.py @@ -0,0 +1,27 @@ +#!/usr/bin/env python + +from lsst.the.monster import ( + GaiaXPMinusDESOffsetMapMaker, + PS1MinusDESOffsetMapMaker, + SkyMapperMinusDESOffsetMapMaker, + VSTMinusDESOffsetMapMaker, + PS1MinusGaiaXPOffsetMapMaker, + SkyMapperMinusGaiaXPOffsetMapMaker, + SkyMapperMinusPS1OffsetMapMaker, + VSTMinusGaiaXPOffsetMapMaker, +) + + +for maker_class in ( + GaiaXPMinusDESOffsetMapMaker, + PS1MinusDESOffsetMapMaker, + SkyMapperMinusDESOffsetMapMaker, + VSTMinusDESOffsetMapMaker, + PS1MinusGaiaXPOffsetMapMaker, + SkyMapperMinusGaiaXPOffsetMapMaker, + SkyMapperMinusPS1OffsetMapMaker, + VSTMinusGaiaXPOffsetMapMaker, +): + maker = maker_class() + hsp_file = maker.measure_offset_map() + maker.plot_offset_maps(hsp_file) diff --git a/python/lsst/the/monster/__init__.py b/python/lsst/the/monster/__init__.py index d2ff64b..26c5c25 100644 --- a/python/lsst/the/monster/__init__.py +++ b/python/lsst/the/monster/__init__.py @@ -25,4 +25,5 @@ from .splinecolorterms import * from .refcats import * from .measure_colorterms import * +from .measure_offsetmaps import * from .utils import * diff --git a/python/lsst/the/monster/measure_offsetmaps.py b/python/lsst/the/monster/measure_offsetmaps.py new file mode 100644 index 0000000..4050dd6 --- /dev/null +++ b/python/lsst/the/monster/measure_offsetmaps.py @@ -0,0 +1,321 @@ +import numpy as np +from astropy import units +import esutil +import hpgeom as hpg +import healsparse as hsp +import skyproj +import scipy.optimize +import matplotlib.pyplot as plt + +import lsst.sphgeom as sphgeom + +from .utils import read_stars +from .refcats import GaiaXPInfo, DESInfo, SkyMapperInfo, PS1Info, VSTInfo + + +__all__ = [ + "OffsetMapMaker", + "GaiaXPMinusDESOffsetMapMaker", + "PS1MinusDESOffsetMapMaker", + "SkyMapperMinusDESOffsetMapMaker", + "VSTMinusDESOffsetMapMaker", + "PS1MinusGaiaXPOffsetMapMaker", + "SkyMapperMinusGaiaXPOffsetMapMaker", + "SkyMapperMinusPS1OffsetMapMaker", + "VSTMinusGaiaXPOffsetMapMaker", +] + + +class OffsetMapMaker: + MinuendInfoClass = None + SubtrahendInfoClass = None + + # Does this have low Galactic latitude area? + has_low_glat = False + + # Name of ID to use for matching. + match_id_name = "GaiaDR3_id" + + @property + def nside(self): + return 128 + + @property + def nside_coarse(self): + return 8 + + @property + def htm_level(self): + return 7 + + def measure_offset_map(self, bands=None, overwrite=False): + """Measure the offset map, and save to a healsparse file. + + Parameters + ---------- + bands : `list` [`str`] + Name of bands to compute. If not specified, will use + the overlap in the two catalogs. + overwrite : `bool`, optional + Overwrite an existing map file. + + Returns + ------- + hsp_file : `str` + Name of healsparse file with maps. + """ + minuend_info = self.MinuendInfoClass() + subtrahend_info = self.SubtrahendInfoClass() + + minuend_path = minuend_info.write_path + subtrahend_path = subtrahend_info.write_path + + if bands is None: + minuend_bands = set(minuend_info.bands) + subtrahend_bands = set(subtrahend_info.bands) + bands = sorted(list(minuend_bands.intersection(subtrahend_bands))) + + print("Using bands: ", bands) + + southern_pixels = hpg.query_box(self.nside_coarse, 0.0, 360.0, -90.0, 35.0) + + healpix_pixelization = sphgeom.HealpixPixelization(hpg.nside_to_order(self.nside_coarse)) + htm_pixelization = sphgeom.HtmPixelization(self.htm_level) + + dtype = [("nmatch", "i4"),] + for band in bands: + dtype.extend([ + (f"offset_{band}", "f4"), + (f"ngood_{band}", "i4"), + ]) + + offset_map = hsp.HealSparseMap.make_empty(32, self.nside, dtype, primary="nmatch") + + for southern_pixel in southern_pixels: + healpix_poly = healpix_pixelization.pixel(southern_pixel) + + htm_pixel_range = htm_pixelization.envelope(healpix_poly) + htm_pixel_list = [] + for r in htm_pixel_range.ranges(): + htm_pixel_list.extend(range(r[0], r[1])) + + minuend_stars = read_stars(minuend_path, htm_pixel_list, allow_missing=True) + if len(minuend_stars) == 0: + continue + subtrahend_stars = read_stars(subtrahend_path, htm_pixel_list, allow_missing=True) + if len(subtrahend_stars) == 0: + continue + + print("Working on ", southern_pixel) + a, b = esutil.numpy_util.match( + minuend_stars[self.match_id_name], + subtrahend_stars[self.match_id_name], + ) + + if len(a) == 0: + continue + + minuend_stars = minuend_stars[a] + subtrahend_stars = subtrahend_stars[b] + + # Cut down to those in the healpix pixel + use, = np.where(hpg.angle_to_pixel( + self.nside_coarse, + subtrahend_stars["coord_ra"], + subtrahend_stars["coord_dec"]) == southern_pixel) + if use.size == 0: + continue + minuend_stars = minuend_stars[use] + subtrahend_stars = subtrahend_stars[use] + + ipnest = hpg.angle_to_pixel( + self.nside, + subtrahend_stars["coord_ra"], + subtrahend_stars["coord_dec"], + ) + + h, rev = esutil.stat.histogram(ipnest, rev=True) + + pixuse, = np.where(h > 0) + + for pixind in pixuse: + i1a = rev[rev[pixind]: rev[pixind + 1]] + + element = np.zeros(1, dtype=dtype) + element["nmatch"] = len(i1a) + + for band in bands: + minuend_flux = minuend_stars[f"decam_{band}_from_{minuend_info.name}_flux"][i1a] + minuend_flux_err = minuend_stars[f"decam_{band}_from_{minuend_info.name}_fluxErr"][i1a] + subtrahend_flux = subtrahend_stars[f"decam_{band}_from_{subtrahend_info.name}_flux"][i1a] + subtrahend_flux_err = subtrahend_stars[ + f"decam_{band}_from_{subtrahend_info.name}_fluxErr" + ][i1a] + + gd = (np.isfinite(minuend_flux) + & np.isfinite(minuend_flux_err) + & np.isfinite(subtrahend_flux) + & np.isfinite(subtrahend_flux_err)) + + element[f"ngood_{band}"] = gd.sum() + + if gd.sum() == 0: + continue + + element[f"offset_{band}"] = np.median( + minuend_flux[gd].quantity.to_value(units.ABmag) + - subtrahend_flux[gd].quantity.to_value(units.ABmag) + ) + + offset_map[ipnest[i1a[0]]] = element + + fname = f"offset_map_{minuend_info.name}-{subtrahend_info.name}.hsp" + offset_map.write(fname) + + return fname + + def plot_offset_maps(self, hsp_file): + """Plot offset maps (and histograms) from a map file. + + Parameters + ---------- + hsp_file : `str` + Name of healsparse file with maps. + """ + minuend_info = self.MinuendInfoClass() + subtrahend_info = self.SubtrahendInfoClass() + + offset_map = hsp.HealSparseMap.read(hsp_file) + + bands = [] + for name in offset_map.dtype.names: + if name.startswith("ngood"): + parts = name.split("_") + bands.append(parts[1]) + + def gauss(x, *p): + A, mu, sigma = p + return A*np.exp(-(x-mu)**2./(2.*sigma**2.)) + + for band in bands: + ngood_band = offset_map[f"ngood_{band}"].copy() + offset_band = offset_map[f"offset_{band}"].copy() + + valid_pixels = ngood_band.valid_pixels + bad, = np.where(ngood_band[valid_pixels] < 3) + offset_band[valid_pixels[bad]] = None + + valid_pixels, valid_ra, valid_dec = offset_band.valid_pixels_pos(return_pixels=True) + + if self.has_low_glat: + # Make a full sky map + + plt.clf() + fig = plt.figure(figsize=(18, 6)) + ax = fig.add_subplot(111) + + sp = skyproj.McBrydeSkyproj(ax=ax) + sp.draw_hspmap(offset_band*1000., zoom=True) + sp.draw_colorbar(label=f"{minuend_info.name} - {subtrahend_info.name} {band} (mmag)") + fig.savefig(f"{minuend_info.name}-{subtrahend_info.name}_fullmap_{band}.png") + plt.close(fig) + + # Now cut out low Galactic latitude regions. + l, b = esutil.coords.eq2gal(valid_ra, valid_dec) + low, = np.where(np.abs(b) < 30.0) + + offset_band[valid_pixels[low]] = None + + plt.clf() + fig = plt.figure(figsize=(18, 6)) + ax = fig.add_subplot(111) + + sp = skyproj.McBrydeSkyproj(ax=ax) + sp.draw_hspmap(offset_band*1000., zoom=True) + sp.draw_colorbar(label=f"{minuend_info.name} - {subtrahend_info.name} {band} (mmag)") + fig.savefig(f"{minuend_info.name}-{subtrahend_info.name}_highglat_{band}.png") + plt.close(fig) + + # Do a histogram and fit a Gaussian to it. + plt.clf() + data = offset_band[offset_band.valid_pixels]*1000. + vmin, vmax = np.percentile(data, q=[1.0, 99.0]) + nbins = 100 + n, b, p = plt.hist( + data, + bins=np.linspace(vmin, vmax, nbins), + histtype="step", + color="blue", + lw=1.5, + ) + + p0 = [data.size, np.mean(data), np.std(data)] + hist_fit_x = (np.array(b[0: -1]) + np.array(b[1:]))/2. + hist_fit_y = np.array(n) + try: + coeff, var_mat = scipy.optimize.curve_fit(gauss, hist_fit_x, hist_fit_y, p0=p0) + except RuntimeError: + # During testing, this may not get a fit. + continue + + xvals = np.linspace(vmin, vmax, 1000) + yvals = gauss(xvals, *coeff) + + plt.plot(xvals, yvals, "k--", linewidth=3) + plt.xlabel(f"{minuend_info.name} - {subtrahend_info.name} {band} (mmag)") + plt.ylabel(f"Number of nside={offset_band.nside_sparse} pixels") + plt.annotate( + r"$\sigma = %.1f\,\mathrm{mmag}$" % (coeff[2]), + (0.99, 0.99), + fontsize=14, + xycoords="axes fraction", + ha="right", + va="top", + ) + plt.savefig(f"{minuend_info.name}-{subtrahend_info.name}_highglat_hist_{band}.png") + + +class GaiaXPMinusDESOffsetMapMaker(OffsetMapMaker): + MinuendInfoClass = GaiaXPInfo + SubtrahendInfoClass = DESInfo + + +class PS1MinusDESOffsetMapMaker(OffsetMapMaker): + MinuendInfoClass = PS1Info + SubtrahendInfoClass = DESInfo + + +class SkyMapperMinusDESOffsetMapMaker(OffsetMapMaker): + MinuendInfoClass = SkyMapperInfo + SubtrahendInfoClass = DESInfo + + +class VSTMinusDESOffsetMapMaker(OffsetMapMaker): + MinuendInfoClass = VSTInfo + SubtrahendInfoClass = DESInfo + + +class PS1MinusGaiaXPOffsetMapMaker(OffsetMapMaker): + MinuendInfoClass = PS1Info + SubtrahendInfoClass = GaiaXPInfo + + has_low_glat = True + + +class SkyMapperMinusGaiaXPOffsetMapMaker(OffsetMapMaker): + MinuendInfoClass = SkyMapperInfo + SubtrahendInfoClass = GaiaXPInfo + + has_low_glat = True + + +class SkyMapperMinusPS1OffsetMapMaker(OffsetMapMaker): + MinuendInfoClass = SkyMapperInfo + SubtrahendInfoClass = PS1Info + + has_low_glat = True + + +class VSTMinusGaiaXPOffsetMapMaker(OffsetMapMaker): + MinuendInfoClass = VSTInfo + SubtrahendInfoClass = GaiaXPInfo From 49541c56a7fbde5c3c8508fb23db1c25cca05700 Mon Sep 17 00:00:00 2001 From: Eli Rykoff Date: Thu, 19 Oct 2023 19:15:53 -0700 Subject: [PATCH 05/11] Add tests to make and plot offset maps. --- tests/test_make_maps.py | 68 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 68 insertions(+) create mode 100644 tests/test_make_maps.py diff --git a/tests/test_make_maps.py b/tests/test_make_maps.py new file mode 100644 index 0000000..0c43f61 --- /dev/null +++ b/tests/test_make_maps.py @@ -0,0 +1,68 @@ +import unittest +import os +import tempfile + +# Ensure that matplotlib doesn't try to open a display during testing. +import matplotlib +matplotlib.use("Agg") + +import lsst.utils # noqa: E402 + + +from lsst.the.monster import PS1MinusGaiaXPOffsetMapMaker # noqa: E402 +from lsst.the.monster import MatchAndTransform # noqa: E402 +from test_catalog_measurement import GaiaDR3InfoTester, GaiaXPInfoTester, PS1InfoTester # noqa: E402 + +ROOT = os.path.abspath(os.path.dirname(__file__)) + + +class MonsterMakeOffsetMapsTest(lsst.utils.tests.TestCase): + def test_MakeOffsetMapAndPlot(self): + # Set up the test. + htmid = 147267 + + gaia_xp_info = GaiaXPInfoTester() + ps1_info = PS1InfoTester() + + with tempfile.TemporaryDirectory() as temp_dir: + # Use yet another way of hacking in test paths + PS1InfoTester.WRITE_PATH = os.path.join(temp_dir, f"{gaia_xp_info.name}_transformed") + GaiaXPInfoTester.WRITE_PATH = os.path.join(temp_dir, f"{ps1_info.name}_transformed") + + class PS1MinusGaiaXPOffsetMapMakerTester(PS1MinusGaiaXPOffsetMapMaker): + MinuendInfoClass = PS1InfoTester + SubtrahendInfoClass = GaiaXPInfoTester + + match_id_name = "TestGaiaDR3_id" + + os.chdir(temp_dir) + + # We first need to do the match and transform. + mat = MatchAndTransform( + catalog_info_class_list=[GaiaXPInfoTester, PS1InfoTester], + gaia_reference_class=GaiaDR3InfoTester, + ) + mat.run(htmid=htmid) + + maker = PS1MinusGaiaXPOffsetMapMakerTester() + hsp_file = maker.measure_offset_map() + + self.assertTrue(os.path.isfile(hsp_file)) + + maker.plot_offset_maps(hsp_file) + + for band in ["g", "r", "i", "z", "y"]: + self.assertTrue(os.path.isfile(f"TestPS1-TestGaiaXP_fullmap_{band}.png")) + self.assertTrue(os.path.isfile(f"TestPS1-TestGaiaXP_highglat_{band}.png")) + # This one fails the fit, which is fine, it's a terrible fit. + if band != "z": + self.assertTrue(os.path.isfile(f"TestPS1-TestGaiaXP_highglat_hist_{band}.png")) + + +def setup_module(module): + lsst.utils.tests.init() + + +if __name__ == "__main__": + lsst.utils.tests.init() + unittest.main() From bee66efb082005e705f6c1d38f36c9f19e8b98c7 Mon Sep 17 00:00:00 2001 From: Eli Rykoff Date: Thu, 19 Oct 2023 19:39:05 -0700 Subject: [PATCH 06/11] Add bands to DESInfo. --- python/lsst/the/monster/refcats.py | 1 + 1 file changed, 1 insertion(+) diff --git a/python/lsst/the/monster/refcats.py b/python/lsst/the/monster/refcats.py index 7183717..e25d945 100644 --- a/python/lsst/the/monster/refcats.py +++ b/python/lsst/the/monster/refcats.py @@ -304,6 +304,7 @@ def get_mag_range(self, band): class DESInfo(RefcatInfo): PATH = "/sdf/data/rubin/shared/the_monster/sharded_refcats/des_y6_calibration_stars_20230511" NAME = "DES" + bands = ["g", "r", "i", "z", "y"] def get_flux_field(self, band): return f"MAG_STD_{band.upper()}_flux" From 32daf22c54609d503b914fda65533230e1b5556d Mon Sep 17 00:00:00 2001 From: Eli Rykoff Date: Fri, 20 Oct 2023 09:23:00 -0700 Subject: [PATCH 07/11] Add additional offset map logging. --- python/lsst/the/monster/measure_offsetmaps.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/python/lsst/the/monster/measure_offsetmaps.py b/python/lsst/the/monster/measure_offsetmaps.py index 4050dd6..463a041 100644 --- a/python/lsst/the/monster/measure_offsetmaps.py +++ b/python/lsst/the/monster/measure_offsetmaps.py @@ -67,6 +67,8 @@ def measure_offset_map(self, bands=None, overwrite=False): minuend_info = self.MinuendInfoClass() subtrahend_info = self.SubtrahendInfoClass() + print(f"Computing offset maps for {minuend_info.name} - {subtrahend_info.name}") + minuend_path = minuend_info.write_path subtrahend_path = subtrahend_info.write_path From 781972d3659d2a6ac1d20e6803866899fb5c8e84 Mon Sep 17 00:00:00 2001 From: Eli Rykoff Date: Fri, 20 Oct 2023 14:56:43 -0700 Subject: [PATCH 08/11] Add ability to remake plots, and ra/dec limits. --- python/lsst/the/monster/measure_offsetmaps.py | 55 +++++++++++++++++-- 1 file changed, 51 insertions(+), 4 deletions(-) diff --git a/python/lsst/the/monster/measure_offsetmaps.py b/python/lsst/the/monster/measure_offsetmaps.py index 463a041..a3c94b6 100644 --- a/python/lsst/the/monster/measure_offsetmaps.py +++ b/python/lsst/the/monster/measure_offsetmaps.py @@ -1,3 +1,4 @@ +import os import numpy as np from astropy import units import esutil @@ -36,6 +37,11 @@ class OffsetMapMaker: # Name of ID to use for matching. match_id_name = "GaiaDR3_id" + # When plotting, limit to pixels within this box. + # (The ra range can wrap around 0.0) + plot_ra_range = [0.0, 360.0] + plot_dec_range = [-90.0, 90.0] + @property def nside(self): return 128 @@ -67,6 +73,15 @@ def measure_offset_map(self, bands=None, overwrite=False): minuend_info = self.MinuendInfoClass() subtrahend_info = self.SubtrahendInfoClass() + fname = f"offset_map_{minuend_info.name}-{subtrahend_info.name}.hsp" + + if os.path.isfile(fname): + if overwrite: + print(f"Found existing {fname}; will overwrite.") + else: + print(f"Found existing {fname}; overwrite=False so nothing to do.") + return fname + print(f"Computing offset maps for {minuend_info.name} - {subtrahend_info.name}") minuend_path = minuend_info.write_path @@ -171,8 +186,7 @@ def measure_offset_map(self, bands=None, overwrite=False): offset_map[ipnest[i1a[0]]] = element - fname = f"offset_map_{minuend_info.name}-{subtrahend_info.name}.hsp" - offset_map.write(fname) + offset_map.write(fname, clobber=overwrite) return fname @@ -187,10 +201,10 @@ def plot_offset_maps(self, hsp_file): minuend_info = self.MinuendInfoClass() subtrahend_info = self.SubtrahendInfoClass() - offset_map = hsp.HealSparseMap.read(hsp_file) + offset_map_in = hsp.HealSparseMap.read(hsp_file) bands = [] - for name in offset_map.dtype.names: + for name in offset_map_in.dtype.names: if name.startswith("ngood"): parts = name.split("_") bands.append(parts[1]) @@ -199,6 +213,18 @@ def gauss(x, *p): A, mu, sigma = p return A*np.exp(-(x-mu)**2./(2.*sigma**2.)) + # Exclude pixels that are outside of the plotting range. + good_pixels = hpg.query_box( + offset_map_in.nside_sparse, + self.plot_ra_range[0], + self.plot_ra_range[1], + self.plot_dec_range[0], + self.plot_dec_range[1], + ) + + offset_map = hsp.HealSparseMap.make_empty_like(offset_map_in) + offset_map[good_pixels] = offset_map_in[good_pixels] + for band in bands: ngood_band = offset_map[f"ngood_{band}"].copy() offset_band = offset_map[f"offset_{band}"].copy() @@ -281,21 +307,38 @@ class GaiaXPMinusDESOffsetMapMaker(OffsetMapMaker): MinuendInfoClass = GaiaXPInfo SubtrahendInfoClass = DESInfo + # Make sure we exclude the DES outrigger fields. + plot_ra_range = [270.0, 120.0] + plot_dec_range = [-80.0, 10.0] + class PS1MinusDESOffsetMapMaker(OffsetMapMaker): MinuendInfoClass = PS1Info SubtrahendInfoClass = DESInfo + # Make sure we exclude the DES outrigger fields + # and anything below dec < -30.0 for PS1. + plot_ra_range = [270.0, 120.0] + plot_dec_range = [-30.0, 10.0] + class SkyMapperMinusDESOffsetMapMaker(OffsetMapMaker): MinuendInfoClass = SkyMapperInfo SubtrahendInfoClass = DESInfo + # Make sure we exclude the DES outrigger fields. + plot_ra_range = [270.0, 120.0] + plot_dec_range = [-80.0, 10.0] + class VSTMinusDESOffsetMapMaker(OffsetMapMaker): MinuendInfoClass = VSTInfo SubtrahendInfoClass = DESInfo + # Make sure we exclude the DES outrigger fields. + plot_ra_range = [270.0, 120.0] + plot_dec_range = [-80.0, 10.0] + class PS1MinusGaiaXPOffsetMapMaker(OffsetMapMaker): MinuendInfoClass = PS1Info @@ -303,6 +346,10 @@ class PS1MinusGaiaXPOffsetMapMaker(OffsetMapMaker): has_low_glat = True + # Make sure we exclude + # anything below dec < -30.0 for PS1. + plot_dec_range = [-30.0, 90.0] + class SkyMapperMinusGaiaXPOffsetMapMaker(OffsetMapMaker): MinuendInfoClass = SkyMapperInfo From 6e72b1f311655c1e9b7dd29d5c007939290e4e59 Mon Sep 17 00:00:00 2001 From: Eli Rykoff Date: Fri, 20 Oct 2023 14:58:37 -0700 Subject: [PATCH 09/11] Improve logging. --- python/lsst/the/monster/measure_offsetmaps.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/python/lsst/the/monster/measure_offsetmaps.py b/python/lsst/the/monster/measure_offsetmaps.py index a3c94b6..217c986 100644 --- a/python/lsst/the/monster/measure_offsetmaps.py +++ b/python/lsst/the/monster/measure_offsetmaps.py @@ -79,7 +79,7 @@ def measure_offset_map(self, bands=None, overwrite=False): if overwrite: print(f"Found existing {fname}; will overwrite.") else: - print(f"Found existing {fname}; overwrite=False so nothing to do.") + print(f"Found existing {fname}; overwrite=False so no need to remake map.") return fname print(f"Computing offset maps for {minuend_info.name} - {subtrahend_info.name}") @@ -201,6 +201,8 @@ def plot_offset_maps(self, hsp_file): minuend_info = self.MinuendInfoClass() subtrahend_info = self.SubtrahendInfoClass() + print(f"Plotting offset maps for {minuend_info.name} - {subtrahend_info.name}") + offset_map_in = hsp.HealSparseMap.read(hsp_file) bands = [] From b7047f7fc4733f0bce6eb3bc171b7fe31c327909 Mon Sep 17 00:00:00 2001 From: Eli Rykoff Date: Fri, 20 Oct 2023 15:11:31 -0700 Subject: [PATCH 10/11] Update some plot settings. --- python/lsst/the/monster/measure_offsetmaps.py | 13 +++++++++++-- 1 file changed, 11 insertions(+), 2 deletions(-) diff --git a/python/lsst/the/monster/measure_offsetmaps.py b/python/lsst/the/monster/measure_offsetmaps.py index 217c986..84a838e 100644 --- a/python/lsst/the/monster/measure_offsetmaps.py +++ b/python/lsst/the/monster/measure_offsetmaps.py @@ -41,6 +41,7 @@ class OffsetMapMaker: # (The ra range can wrap around 0.0) plot_ra_range = [0.0, 360.0] plot_dec_range = [-90.0, 90.0] + plot_figsize = (16, 6) @property def nside(self): @@ -241,7 +242,7 @@ def gauss(x, *p): # Make a full sky map plt.clf() - fig = plt.figure(figsize=(18, 6)) + fig = plt.figure(figsize=self.plot_figsize) ax = fig.add_subplot(111) sp = skyproj.McBrydeSkyproj(ax=ax) @@ -257,7 +258,7 @@ def gauss(x, *p): offset_band[valid_pixels[low]] = None plt.clf() - fig = plt.figure(figsize=(18, 6)) + fig = plt.figure(figsize=self.plot_figsize) ax = fig.add_subplot(111) sp = skyproj.McBrydeSkyproj(ax=ax) @@ -312,6 +313,7 @@ class GaiaXPMinusDESOffsetMapMaker(OffsetMapMaker): # Make sure we exclude the DES outrigger fields. plot_ra_range = [270.0, 120.0] plot_dec_range = [-80.0, 10.0] + plot_figsize = (10, 6) class PS1MinusDESOffsetMapMaker(OffsetMapMaker): @@ -322,6 +324,7 @@ class PS1MinusDESOffsetMapMaker(OffsetMapMaker): # and anything below dec < -30.0 for PS1. plot_ra_range = [270.0, 120.0] plot_dec_range = [-30.0, 10.0] + plot_figsize = (12, 6) class SkyMapperMinusDESOffsetMapMaker(OffsetMapMaker): @@ -331,6 +334,7 @@ class SkyMapperMinusDESOffsetMapMaker(OffsetMapMaker): # Make sure we exclude the DES outrigger fields. plot_ra_range = [270.0, 120.0] plot_dec_range = [-80.0, 10.0] + plot_figsize = (10, 6) class VSTMinusDESOffsetMapMaker(OffsetMapMaker): @@ -340,6 +344,7 @@ class VSTMinusDESOffsetMapMaker(OffsetMapMaker): # Make sure we exclude the DES outrigger fields. plot_ra_range = [270.0, 120.0] plot_dec_range = [-80.0, 10.0] + plot_figsize = (12, 6) class PS1MinusGaiaXPOffsetMapMaker(OffsetMapMaker): @@ -366,6 +371,10 @@ class SkyMapperMinusPS1OffsetMapMaker(OffsetMapMaker): has_low_glat = True + # Make sure we exclude + # anything below dec < -30.0 for PS1. + plot_dec_range = [-30.0, 90.0] + class VSTMinusGaiaXPOffsetMapMaker(OffsetMapMaker): MinuendInfoClass = VSTInfo From e0b95181f2a6490bce6534ddf9b428475f6b1293 Mon Sep 17 00:00:00 2001 From: Eli Rykoff Date: Fri, 20 Oct 2023 15:38:01 -0700 Subject: [PATCH 11/11] Add offset to annotation. --- python/lsst/the/monster/measure_offsetmaps.py | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/python/lsst/the/monster/measure_offsetmaps.py b/python/lsst/the/monster/measure_offsetmaps.py index 84a838e..7506328 100644 --- a/python/lsst/the/monster/measure_offsetmaps.py +++ b/python/lsst/the/monster/measure_offsetmaps.py @@ -296,13 +296,22 @@ def gauss(x, *p): plt.xlabel(f"{minuend_info.name} - {subtrahend_info.name} {band} (mmag)") plt.ylabel(f"Number of nside={offset_band.nside_sparse} pixels") plt.annotate( - r"$\sigma = %.1f\,\mathrm{mmag}$" % (coeff[2]), + r"$\mu = %.1f\,\mathrm{mmag}$" % (coeff[1]), (0.99, 0.99), fontsize=14, xycoords="axes fraction", ha="right", va="top", ) + plt.annotate( + r"$\sigma = %.1f\,\mathrm{mmag}$" % (coeff[2]), + (0.99, 0.92), + fontsize=14, + xycoords="axes fraction", + ha="right", + va="top", + ) + plt.savefig(f"{minuend_info.name}-{subtrahend_info.name}_highglat_hist_{band}.png")