Skip to content

Aperture photometry #371

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 7 commits into
base: photometry-dev
Choose a base branch
from
Open

Aperture photometry #371

Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
e6a058f
Minimal updates to get notebook working
mwcraig May 5, 2024
b5227c1
Move last STScI notebook into main photometry docs
mwcraig May 5, 2024
731c25a
Minor cosmetic fixes
mwcraig Jul 15, 2024
22064ab
Use angular scale for sky apertures
mwcraig Jul 15, 2024
9658484
Use local file instead of URL
mwcraig Jul 15, 2024
5604a79
Use high level WCS API
mwcraig Jul 15, 2024
8a39dd3
Format colorbar properly
mwcraig Jul 15, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
89 changes: 50 additions & 39 deletions ...e_photometry/03_aperture_photometry.ipynb → ...ry/04.01.02-aperture-photometry-xdf.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -127,7 +127,8 @@
"outputs": [],
"source": [
"url = 'https://archive.stsci.edu/pub/hlsp/xdf/hlsp_xdf_hst_acswfc-60mas_hudf_f435w_v1_sci.fits'\n",
"xdf_image = CCDData.read(url, unit=u.electron / u.s)"
"file = 'hlsp_xdf_hst_acswfc-60mas_hudf_f435w_v1_sci.fits'\n",
"xdf_image = CCDData.read(file, unit=u.electron / u.s)"
]
},
{
Expand Down Expand Up @@ -158,9 +159,7 @@
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"scrolled": false
},
"metadata": {},
"outputs": [],
"source": [
"# Set up the figure with subplots\n",
Expand Down Expand Up @@ -226,7 +225,7 @@
"metadata": {},
"outputs": [],
"source": [
"from photutils import find_peaks\n",
"from photutils.detection import find_peaks\n",
"from photutils.centroids import centroid_2dg"
]
},
Expand All @@ -247,7 +246,7 @@
"outputs": [],
"source": [
"sources_findpeaks = find_peaks(xdf_image.data, mask=xdf_image.mask, \n",
" threshold=20.*std, box_size=30, \n",
" threshold=20.*std, box_size=29, \n",
" centroid_func=centroid_2dg) \n",
"# Display the table\n",
"sources_findpeaks"
Expand Down Expand Up @@ -302,7 +301,7 @@
"metadata": {},
"outputs": [],
"source": [
"from photutils import CircularAperture"
"from photutils.aperture import CircularAperture, EllipticalAperture, EllipticalAnnulus"
]
},
{
Expand All @@ -312,7 +311,7 @@
"outputs": [],
"source": [
"# Define the aperture\n",
"position = (sources_findpeaks['x_centroid'], sources_findpeaks['y_centroid'])\n",
"position = list((x, y) for x, y in zip(sources_findpeaks['x_centroid'], sources_findpeaks['y_centroid']))\n",
"radius = 10. # pixels\n",
"circular_aperture = CircularAperture(position, r=radius)"
]
Expand Down Expand Up @@ -371,8 +370,10 @@
"metadata": {},
"outputs": [],
"source": [
"from photutils import (detect_sources, source_properties, \\\n",
" EllipticalAnnulus, EllipticalAperture)"
"# from photutils import (detect_sources, source_properties, \\\n",
"# EllipticalAnnulus, EllipticalAperture)\n",
"\n",
"from photutils.segmentation import detect_sources, SourceCatalog"
]
},
{
Expand All @@ -387,9 +388,7 @@
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"scrolled": false
},
"metadata": {},
"outputs": [],
"source": [
"# Define threshold and minimum object size\n",
Expand All @@ -400,7 +399,7 @@
"segm = detect_sources(xdf_image.data, threshold, npixels)\n",
"\n",
"# Create a catalog using source properties\n",
"catalog = source_properties(xdf_image.data, segm)\n",
"catalog = SourceCatalog(xdf_image.data, segm)\n",
"table = catalog.to_table()\n",
"\n",
"# Display the table\n",
Expand Down Expand Up @@ -437,9 +436,9 @@
"# Create the apertures\n",
"elliptical_apertures = []\n",
"for obj in catalog:\n",
" position = (obj.xcentroid.value, obj.ycentroid.value)\n",
" a = obj.semimajor_axis_sigma.value * r # pixel\n",
" b = obj.semiminor_axis_sigma.value * r # pixel\n",
" position = (obj.xcentroid, obj.ycentroid)\n",
" a = obj.semimajor_sigma.value * r # pixel\n",
" b = obj.semiminor_sigma.value * r # pixel\n",
" theta = obj.orientation.to(u.rad).value\n",
" elliptical_apertures.append(EllipticalAperture(position, a, b, theta=theta))"
]
Expand Down Expand Up @@ -521,7 +520,7 @@
"metadata": {},
"outputs": [],
"source": [
"from photutils import SkyEllipticalAperture\n",
"from photutils.aperture import SkyEllipticalAperture\n",
"from astropy.coordinates import SkyCoord\n",
"\n",
"r = 3. # pixels; approximate isophotal extent of semimajor axis\n",
Expand All @@ -530,13 +529,18 @@
"sky_elliptical_apertures = []\n",
"for obj in catalog:\n",
" # Convert the centroids into RA/Dec using WCS\n",
" ra, dec = wcs.all_pix2world(obj.xcentroid.value, obj.ycentroid.value, 0)\n",
" coord = wcs.pixel_to_world_values([[obj.xcentroid, obj.ycentroid]])\n",
" ra, dec = coord[0]\n",
" # Convert the positions to an Astropy SkyCoord object, with units!\n",
" sky_position = SkyCoord(ra, dec, unit=u.deg)\n",
" \n",
"\n",
" # NOT 100% SURE THIS IS CORRECT -- PROJ PLANE IS NOT IMAGE PLANE\n",
" # Get the pixel scale, i.e. the angular size of a pixel, from the wcs.\n",
" pixel_scale_x, pixel_scale_y = wcs.proj_plane_pixel_scales()\n",
"\n",
" # Define the elliptical parameters, now with units\n",
" a = obj.semimajor_axis_sigma.value * r * u.pix\n",
" b = obj.semiminor_axis_sigma.value * r * u.pix\n",
" a = obj.semimajor_sigma.value * r * pixel_scale_x\n",
" b = obj.semiminor_sigma.value * r * pixel_scale_y\n",
" theta = obj.orientation.value * u.rad\n",
" \n",
" # Convert the theta from radians from X axis to the radians from North \n",
Expand Down Expand Up @@ -592,7 +596,7 @@
"metadata": {},
"outputs": [],
"source": [
"from photutils import aperture_photometry\n",
"from photutils.aperture import aperture_photometry\n",
"from astropy.table import QTable"
]
},
Expand Down Expand Up @@ -646,7 +650,7 @@
"\n",
"Since our elliptical apertures are distinct apertures at distinct positions, we need to do a little more work to get a single table of photometric values.\n",
"\n",
"(This step will take a while, almost 5 minutes, so hang tight!)"
"(This step will take a while, almost NO IT DOESN'T 5 minutes, so hang tight!)"
]
},
{
Expand Down Expand Up @@ -711,9 +715,7 @@
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"scrolled": false
},
"metadata": {},
"outputs": [],
"source": [
"plt.figure(figsize=(8, 5))\n",
Expand All @@ -724,7 +726,7 @@
"plt.xscale('log')\n",
"plt.title('Count Rate v. Aperture Area')\n",
"plt.xlabel('Aperture Area [pixels$^2$]')\n",
"plt.ylabel(r'Flux Count Rate ({})'.format(xdf_image.unit.to_string('latex')))"
"plt.ylabel(r'Flux Count Rate ({})'.format(xdf_image.unit.to_string('latex')));"
]
},
{
Expand Down Expand Up @@ -783,10 +785,10 @@
"# Create the apertures\n",
"elliptical_annuli = []\n",
"for obj in catalog:\n",
" position = (obj.xcentroid.value, obj.ycentroid.value)\n",
" a_in = obj.semimajor_axis_sigma.value * r_in\n",
" a_out = obj.semimajor_axis_sigma.value * r_out\n",
" b_out = obj.semiminor_axis_sigma.value * r_out\n",
" position = (obj.xcentroid, obj.ycentroid)\n",
" a_in = obj.semimajor_sigma.value * r_in\n",
" a_out = obj.semimajor_sigma.value * r_out\n",
" b_out = obj.semiminor_sigma.value * r_out\n",
" theta = obj.orientation.to(u.rad).value\n",
" elliptical_annuli.append(EllipticalAnnulus(position, a_in, a_out, b_out, theta=theta))"
]
Expand All @@ -811,9 +813,18 @@
" aperture.plot(color='white', alpha=0.7, axes=ax1)\n",
"\n",
"# Define the colorbar\n",
"cbar = plt.colorbar(fitsplot, fraction=0.046, pad=0.04, ticks=LogLocator(subs=range(10)))\n",
"labels = ['$10^{-4}$'] + [''] * 8 + ['$10^{-3}$'] + [''] * 8 + ['$10^{-2}$']\n",
"cbar.ax.set_yticklabels(labels)\n",
"cbar = plt.colorbar(fitsplot, fraction=0.046, pad=0.04, ticks=LogLocator())\n",
"\n",
"def format_colorbar(bar):\n",
" # Add minor tickmarks\n",
" bar.ax.yaxis.set_minor_locator(LogLocator(subs=range(1, 10)))\n",
"\n",
" # Force the labels to be displayed as powers of ten and only at exact powers of ten\n",
" bar.ax.set_yticks([1e-4, 1e-3, 1e-2])\n",
" labels = [f'$10^{{{pow:.0f}}}$' for pow in np.log10(bar.ax.get_yticks())]\n",
" bar.ax.set_yticklabels(labels)\n",
"\n",
"format_colorbar(cbar)\n",
"\n",
"# Define labels\n",
"cbar.set_label(r'Flux Count Rate ({})'.format(xdf_image.unit.to_string('latex')), \n",
Expand All @@ -824,7 +835,7 @@
"\n",
"# Crop to show an inset of the data\n",
"ax1.set_xlim(2000, 3000)\n",
"ax1.set_ylim(2000, 1000)"
"ax1.set_ylim(2000, 1000);"
]
},
{
Expand Down Expand Up @@ -987,7 +998,7 @@
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
Expand All @@ -1001,9 +1012,9 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.5"
"version": "3.11.9"
}
},
"nbformat": 4,
"nbformat_minor": 2
"nbformat_minor": 4
}