Run wcs-based matching to compare object tables to metadetect no-shear catalogs
This takes as the metadetect no-shear catalog and object catalogs as inputs and runs matching using the wcs-based WcsMatch matcher.
Standard imports
[1]:
import hpmcm
import tables_io
import glob
import os
import numpy as np
import matplotlib.pyplot as plt
Set up the configuration
[2]:
DATADIR = "test_data" # Input data directory
shear_st = "0p01" # Applied shear as a string
shear = 0.01 # Decimal version of applied shear
shear_type = "wmom" # which object characterization to use
tract = 10463 # which tract to study
REF_DIR = (37.9, 7.0) # RA, DEC in deg of center of match region
REGION_SIZE = (0.375, 0.375) # Size of match region in degrees
PIXEL_SIZE = 0.5/3600. # Size of pixels used in matching
PIXEL_R2CUT = 4. # Cut at distance**2 = 4 pixels
PIXEL_MATCH_SCALE = 1 # Use pixel scale to do matching
SOURCE_TABLEFILES = sorted(glob.glob(os.path.join(DATADIR, f"shear_*_{shear_st}_cleaned_{tract}_ns.pq")))
SOURCE_TABLEFILES.append(os.path.join(DATADIR, f"object_{tract}.pq"))
SOURCE_TABLEFILES.reverse()
SOURCE_TABLEFILES = [SOURCE_TABLEFILES[0], SOURCE_TABLEFILES[1]]
VISIT_IDS = np.arange(len(SOURCE_TABLEFILES))
Make the matcher, reduce the data
[3]:
matcher = hpmcm.WcsMatch.create(REF_DIR, REGION_SIZE, pixel_size=PIXEL_SIZE, pixel_r2_cut=PIXEL_R2CUT)
matcher.reduceData(SOURCE_TABLEFILES, VISIT_IDS)
/home/docs/checkouts/readthedocs.org/user_builds/hpmcm/envs/stable/lib/python3.11/site-packages/hpmcm/wcs_match.py:157: SettingWithCopyWarning:
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead
See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
df_clean["x_pix"] = x_pix
/home/docs/checkouts/readthedocs.org/user_builds/hpmcm/envs/stable/lib/python3.11/site-packages/hpmcm/wcs_match.py:158: SettingWithCopyWarning:
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead
See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
df_clean["y_pix"] = y_pix
/home/docs/checkouts/readthedocs.org/user_builds/hpmcm/envs/stable/lib/python3.11/site-packages/hpmcm/wcs_match.py:157: SettingWithCopyWarning:
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead
See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
df_clean["x_pix"] = x_pix
/home/docs/checkouts/readthedocs.org/user_builds/hpmcm/envs/stable/lib/python3.11/site-packages/hpmcm/wcs_match.py:158: SettingWithCopyWarning:
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead
See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
df_clean["y_pix"] = y_pix
Make a plot comparing the signal-to-noise in the two catalogs
[4]:
_ = plt.hist(matcher.full_data[0].snr, bins=np.logspace(0, 4, 81), alpha=0.4, label="obj")
_ = plt.hist(matcher.full_data[1].snr, bins=np.logspace(0, 4, 81), alpha=0.4, label="pgauss")
_ = plt.xscale('log')
_ = plt.legend()
This should have made 3 x 3 cells
[5]:
matcher.n_cell
[5]:
array([3., 3.])
Run the loop over cells
[6]:
matcher.analysisLoop()
.. Done!
Show a single cluster
The x and y axes here are the in the cluster frame for a single cluster. The color scale shows the number of sources per/pixel.
The x markers are the original source postions. The o makters are the deshear positions.
[7]:
cell = matcher.cell_dict[matcher.getCellIdx(2,2)]
od = cell.analyze(None, 4)
[8]:
cluster = list(cell.cluster_dict.values())[0]
_ = hpmcm.viz_utils.showCluster(od['image'], cluster, cell)
Classify the objects by match type
This looks at the characteristics of the matched objects and categorizes them.
[9]:
obj_lists = hpmcm.classify.classifyObjects(matcher, snr_cut=10)
hpmcm.classify.printObjectTypes(obj_lists)
All Objects: 19023
cut 1 0
cut 2 0
Used: 19023
good (n source from n catalogs): 3710
good faint 3176
faint (< n sources, snr < cut): 11748
mixed (n source from < n catalogs): 12
edge_mixed (mixed near edge of cell): 0
edge_missing (< n sources, near edge of cell): 0
edge_extra (> n sources, near edge of cell): 0
faint (< n sources, snr < cut): 11748
orphan (split off from larger cluster 2
one missing (n-1 sources, not near edge): 350
two missing (n-2 sources, not near edge): 0
many missing (< n-2 sources, not near edge): 0
extra (> n sources, not near edge): 25
[10]:
n_good = len(obj_lists['ideal'])
bad_list = ['edge_mixed', 'edge_missing', 'edge_extra', 'orphan', 'missing', 'two_missing', 'many_missing', 'extra', 'caught']
n_bad = np.sum([len(obj_lists[x]) for x in bad_list])
[11]:
effic = n_good/(n_good+n_bad)
effic_err = np.sqrt(effic*(1-effic)/(n_good+n_bad))
print(f"Effic: {effic:.5} +- {effic_err:.5f}")
Effic: 0.90776 +- 0.00453
Classify objcts using the object table as the reference
[12]:
odict = hpmcm.classify.matchObjectsAgainstRef(matcher, snrCut=10.)
hpmcm.classify.printObjectMatchTypes(odict)
All 19023
Used 19023
New 31
New (faint) 47
In Ref 18945
Faint 11406
Good 4702
Good (faint) 2184
Missing 640
Two Missing 0
All Missing 0
Extra 8
Caught 5
Display an object
[13]:
_ = hpmcm.viz_utils.showShearObj(matcher, obj_lists['missing'][0])
[ ]:
[ ]: