Measuring magnitude at specific RA/Dec points

Hello everyone,

I am currently working on a project involving the injection and extraction of microlensing light curves in DP0. Specifically, I am focusing on measuring the magnitude at specific RA/Dec points across several calexps to recreate a light curve. My goal is to scale this process for a large number of microlensing events.

I have been experimenting with two methods in the LSST framework to achieve this and would greatly appreciate any insights or recommendations.

  1. SingleFrameMeasurementTask: This method provides coherent magnitude and error values consistent with the injected sources. However, it processes all detected sources in the calibrated exposure (calexp), which significantly increases the computation time since I’m only interested in specific points of interest.
    Here are the key highlights from my code that provide context:
# Task configuration
config = SourceDetectionTask.ConfigClass()
config.thresholdValue = 5
config.thresholdType = "stdev"
detect_task = SourceDetectionTask(schema=schema, config=config)
schema = afwTable.SourceTable.makeMinimalSchema()
algMetadata = dafBase.PropertyList()
config = SingleFrameMeasurementTask.ConfigClass()
measure_task = SingleFrameMeasurementTask(schema=schema, config=config,  algMetadata=algMetadata)

# Run tasks
detect_result = detect_task.run(tab, injected_exposure)
sources = detect_result.sources
measure_task.run(measCat=sources, exposure=injected_exposure)

# Finding the points of interest
source = sources[(abs((sources["coord_ra"]/np.pi*180) - rai) < delta * k) & 
                                     (abs((sources["coord_dec"]/np.pi*180) - deci) < delta * k)]
  1. ForcedMeasurementTask: Using this approach, I can measure the flux specifically at my points of interest. However, I am encountering issues with the returned flux values: they are often negative, leading to NaN magnitudes or values greater than expected.
    Here are the key highlights from my code that provide context:
schema = afwTable.SourceTable.makeMinimalSchema()
alias = schema.getAliasMap() 
x_key = schema.addField("centroid_x", type="D")
y_key = schema.addField("centroid_y", type="D")
alias.set("slot_Centroid", "centroid")

xx_key = schema.addField("shape_xx", type="D")
yy_key = schema.addField("shape_yy", type="D")
xy_key = schema.addField("shape_xy", type="D")
alias.set("slot_Shape", "shape")
type_key = schema.addField("type_flag", type="F")
forcedSource = afwTable.SourceCatalog(schema)
config = ForcedMeasurementTask.ConfigClass()
config.copyColumns = {}
config.plugins.names = ['base_SdssCentroid', "base_TransformedCentroid",
    "base_PsfFlux",
    "base_TransformedShape",
]
config.doReplaceWithNoise = False
forcedMeasurementTask = ForcedMeasurementTask(schema, config=config)

config = CharacterizeImageTask.ConfigClass()
config.psfIterations = 1
charImageTask = CharacterizeImageTask(config=config)

inject_config = VisitInjectConfig()
inject_task = VisitInjectTask(config=inject_config)

My questions:

  1. Is it possible to measure magnitudes at specific RA/Dec points without processing all detected sources, using SingleFrameMeasurementTask?
  2. For ForcedMeasurementTask, are there any common pitfalls or parameter settings that might lead to negative flux values? What might be the cause of these NaN or unexpectedly high magnitude readings

Slightly off-topic questions:

  1. Is a RAFT equivalent to a detector in the LSST context? Given the calexp sizes, it seems that RAFT might not equal CCD, but I would appreciate clarification on this.
  2. Regarding the calexp collection using the HTM level: Why is it that some calexps are not captured, even though the HTM level 20 triangles have a boundary circle radius of approximately 0.06 arcsec? Given that a pixel is 0.2 arcsec (which is much larger than 0.06 arcsec), I would expect all points to be covered.
  3. Regarding calexp plotting: Another issue I’m facing: when I plot some calexps, the coordinate axes seem to be offset from the WCS center. Could I be missing something about the astrometry? The approach I used comes from a DP0 tutorial, which relies on astropy and matplotlib.
    Here is the code used for plotting:
# RA (hourangle)
fig = plt.figure()
ax = plt.subplot(projection=WCS(calexp.getWcs().getFitsMetadata()))
calexp_extent = (calexp.getBBox().beginX, calexp.getBBox().endX,
                 calexp.getBBox().beginY, calexp.getBBox().endY)
ax.set_xlabel('RA (hourangle)', fontsize=8)
ax.set_ylabel('Dec (degrees)', fontsize=8)
ax.coords['ra'].set_ticks(number=20)
ax.coords['dec'].set_ticks(number=20)
im = plt.imshow(calexp.image.array, cmap='gray', vmin=-200.0, vmax=400,
                extent=calexp_extent, origin='lower')
plt.grid(color='white', ls='solid')

For example:

  • visit = 211488; detector = 103: This calexp plots correctly centered on the WCS.
  • visit = 214463; detector = 171: In this case, the WCS center is ~60.04 degrees, but the plot shows ~60.18 degrees (shifted by ~0.14 degrees).

Any guidance or tips would be greatly appreciated. Thank you!

1 Like

Hi @nowokaren, thanks for these questions.

Questions #1 and #2

I think it might be possible to pass only a list of source locations to SingeFrameMeasurementTask, but using ForcedMeasurementTask is (I think) the better task for the job of doing photometry at known coordinates. So first let’s focus on your question #2.

Common pitfalls that might lead to NaN or unexpected flux might be when sources are injected near or outside an image’s boundary, or on top of sources in the image. In DP0.2 tutorial notebook 05, Section 5.4 explores this a bit. Here, it would be helpful to explore and characterize the injected sources for which unexpected forced flux measurements are returned. Are they all in a certain region? A certain calexp? Close to other sources? That kind of information would help to diagnose the issue.

Off-topic question #3 I can answer here: a raft is nine CCDs in a 3x3 arrangement.

Off-topic question #4 would be better as its own separate Forum Topic. I don’t quite understand the issue and can’t reproduce it with the information given – but a new Topic would provide more space for, e.g., the code that was not returning the expected results. Potentially relevant here is that use of HTM for spatial queries with the butler has been superseded by the OVERLAPS POINT(ra, dec) functionality demonstrated in, e.g., Section 3.3.2 of DP0.2 tutorial notebook 04b. Maybe check that out and then create a new Topic if the issue in Question #4 persists?

Off-topic question #5 would also be better as it’s own separate Forum Topic. I created “Calexps displayed with imshow seem offset from WCS center” and posted a potential answer there, please respond if it solves the issue?

We can continue to discuss questions #1 and #2 in this thread.

Hi @nowokaren , just wanted to follow up to see if have any responses or further questions on @MelissaGraham’s reply. Sounds like it would be great to continue discussing your questions 1 and 2 here.