Skip to content
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

Jump to bottom

updated with s3df file stream #31

Merged
merged 7 commits into from
Jun 17, 2024
Merged

updated with s3df file stream #31

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
cf9f169
updated with s3df file stream
Jun 12, 2024
16608e8
Merge branch 'main' of github.com:spacetelescope/roman_notebooks into…
Jun 14, 2024
607682e
reordered roman_datamodels and asdf sections
Jun 14, 2024
534c4ad
clear notebook
Jun 15, 2024
aad4415
updated shape-measurement notebook using latest catalog
Jun 15, 2024
d11d843
updated shape-measurement notebook back to S3 bucket files
Jun 17, 2024
8f1bdec
using rdm, fixed some typos/text
Jun 17, 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
113 changes: 71 additions & 42 deletions content/notebooks/measuring_galaxy_shapes/measuring_galaxy_shapes.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -42,7 +42,8 @@
"- *webbpsf* to access Roman's point-spread function.\n",
"- *galsim* to measure source moments and generate Sérsic profiles.\n",
"- *romanisim.bandpass.get_abflux* to convert flux from maggies to photons/s consistently with our simulation.\n",
"- *scipy.optimize.curve_fit* to fit our data."
"- *scipy.optimize.curve_fit* to fit our data.\n",
"- *roman_datamodels* to handle Roman data."
]
},
{
Expand All @@ -63,13 +64,13 @@
"import webbpsf\n",
"import galsim\n",
"import gwcs\n",
"import s3fs\n",
"import copy\n",
"from astropy.coordinates import SkyCoord, Angle\n",
"from astropy.coordinates import SkyCoord\n",
"import astropy.units as u\n",
"from galsim.roman import collecting_area\n",
"from romanisim.bandpass import get_abflux\n",
"import roman_datamodels as rdm\n",
"from scipy.optimize import curve_fit\n",
"import roman_datamodels as rdm\n",
"import s3fs\n",
"%matplotlib inline"
]
},
Expand All @@ -88,7 +89,7 @@
"[REGAUSS method in its HSM module](https://galsim-developers.github.io/GalSim/_build/html/hsm.html) ([Hirata & Seljak 2003](https://ui.adsabs.harvard.edu/abs/2003MNRAS.343..459H/abstract), [Mandelbaum et al. 2005](https://ui.adsabs.harvard.edu/abs/2005MNRAS.361.1287M/abstract)). Second, we fit a galaxy cutout to a Sérsic model.\n",
"\n",
"This notebook relies on Advanced Scientific Data Format (`ASDF`) file manipulations. For additional information about these files\n",
"and how to use them, please check our `exploring_ASDF` notebook.\n",
"and how to use them, please check our `working_with_asdf` notebook.\n",
"\n",
"### Defining terms\n",
"\n",
Expand All @@ -111,7 +112,7 @@
},
"source": [
"## Loading data\n",
"The first step of the analysis is to read the Roman image data, which are stored in ASDF format. For this example we start with a calibrated level-2 simulated image created with [`romanisim`](https://romanisim.readthedocs.io). For more information about Roman's data products check the [Roman User Documentation](roman-docs.stsci.edu)."
"The first step of the analysis is to read the Roman image data, which are stored in ASDF format. For this example we start with a calibrated level-2 simulated image created with [`romanisim`](https://romanisim.readthedocs.io). For more information about Roman's data products check the [Roman User Documentation](roman-docs.stsci.edu). We use `roman_datamodels` to open the simulated image."
]
},
{
Expand All @@ -124,16 +125,13 @@
},
"outputs": [],
"source": [
"# This will change when the data are elsewhere\n",
"asdf_dir_uri = 's3://roman-sci-test-data-prod-summer-beta-test/'\n",
"fs = s3fs.S3FileSystem()\n",
"asdf_file_uri_l2 = asdf_dir_uri + 'ROMANISIM/GALAXIES/r0003201001001001004_01101_0001_WFI01_cal.asdf'\n",
"\n",
"asdf_file_uri = asdf_dir_uri + 'ROMANISIM/GALAXIES/r0003201001001001004_01101_0001_WFI01_cal.asdf'\n",
"\n",
"with fs.open(asdf_file_uri, 'rb') as f:\n",
" dm = rdm.open(f)\n",
" meta = dm.meta.copy()\n",
" image = dm.data.copy()\n",
" wcs = copy.deepcopy(dm.meta.wcs)"
"with fs.open(asdf_file_uri_l2, 'rb') as f:\n",
" img_asdf = rdm.open(f).copy()"
]
},
{
Expand All @@ -151,8 +149,8 @@
"metadata": {},
"outputs": [],
"source": [
"img_arr = image.value\n",
"print('The units of the original data are: ', image.unit)"
"img_arr = img_asdf['data'][:, :].value\n",
"print('The units of the original data are: ', img_asdf['data'].unit)"
]
},
{
Expand All @@ -168,9 +166,9 @@
"metadata": {},
"outputs": [],
"source": [
"catalog_file_uri = asdf_dir_uri + 'ROMANISIM/CATALOGS_SCRIPTS/galfield.ecsv'\n",
"with fs.open(catalog_file_uri, 'rb') as catalog_file_stream:\n",
" catalog = Table.read(catalog_file_stream, format='ascii.ecsv')"
"cat_uri = asdf_dir_uri + 'ROMANISIM/CATALOGS_SCRIPTS/galfield.ecsv'\n",
"with fs.open(cat_uri, 'r') as catalog_file_stream:\n",
" catalog = Table.read(cat_uri, format='ascii.ecsv').copy()"
]
},
{
Expand All @@ -186,12 +184,10 @@
"metadata": {},
"outputs": [],
"source": [
"ra = Angle(catalog['ra'] * u.deg)\n",
"dec = Angle(catalog['dec'] * u.deg)\n",
"ra.wrap_at(180 * u.deg, inplace=True)\n",
"plt.scatter(ra, dec, s=0.1)\n",
"catalog['ra'][catalog['ra'] > 180] -= 360\n",
"plt.scatter(catalog['ra'][::100], catalog['dec'][::100], s=0.1)\n",
"plt.xlabel('RA [deg]', fontsize=16)\n",
"plt.ylabel('Dec [deg]', fontsize=16);"
"plt.ylabel('Dec [deg]', fontsize=16)"
]
},
{
Expand All @@ -215,7 +211,17 @@
"metadata": {},
"outputs": [],
"source": [
"coords = SkyCoord(ra=catalog['ra']*u.deg, dec=catalog['dec']*u.deg) # Save the catalogs' coordinates in a SkyCoord object"
"w = img_asdf['meta']['wcs'] # We call the GWCS object w for later usage"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"coords = SkyCoord(ra=catalog['ra']*u.deg, dec=catalog['dec']*u.deg) # Save the catalogs' coordinates in a SkyCoord object\n",
"y, x = w.world_to_array_index_values(coords) # Convert from sky to detector indices"
]
},
{
Expand All @@ -231,7 +237,11 @@
"metadata": {},
"outputs": [],
"source": [
"gal_sel = (-2.5*np.log10(catalog['F106']) > 21) & (-2.5*np.log10(catalog['F106']) < 21.5) & (catalog['type'] == 'SER')\n",
"gal_sel =(-2.5*np.log10(catalog['F106']) > 21) & (-2.5*np.log10(catalog['F106']) < 21.5) & (catalog['type'] == 'SER')\n",
"# We make sure that the sources are in the sensor\n",
"padding = 10\n",
"inchip = (x > padding) & (x < img_asdf['data'].shape[0] - padding) & (y > padding) & (y < img_asdf['data'].shape[1] - padding)\n",
"gal_sel = (gal_sel) & (inchip)\n",
"print(np.count_nonzero(gal_sel), 'galaxies pass these selection criteria')"
]
},
Expand All @@ -248,7 +258,7 @@
"metadata": {},
"outputs": [],
"source": [
"band = meta.instrument.optical_element\n",
"band = img_asdf['meta']['instrument']['optical_element']\n",
"print('The simulated band is:', band)\n",
"nc = webbpsf.WFI() # We tell webbpsf that we want a Roman PSF\n",
"nc.filter = band # We configure the filter\n",
Expand All @@ -273,7 +283,9 @@
"outputs": [],
"source": [
"# We are going to select one of the galaxies within our galaxy selection gal_sel \n",
"igal = 12\n",
"igal = 0\n",
"ra = catalog['ra'][gal_sel][igal]\n",
"dec = catalog['dec'][gal_sel][igal]\n",
"# We create a fairly large cutout -- 4 times the half-light radius + 1 pixel for padding\n",
"size = int(catalog['half_light_radius'][gal_sel][igal]*4/0.1+1)\n",
"# We cut our original image on the selected coordinates of the source of interest\n",
Expand All @@ -285,7 +297,7 @@
"plt.imshow(cutout.array, origin='lower', norm=LogNorm())\n",
"plt.colorbar(label='DN/s')\n",
"plt.xlabel('X [pix]')\n",
"plt.ylabel('Y [pix]');"
"plt.ylabel('Y [pix]')"
]
},
{
Expand Down Expand Up @@ -377,6 +389,22 @@
"zp = get_abflux(band)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"`galsim` sources interpret their input fluxes as photons/s/cm$^{2}$. So we compute the zeropoint shift from AB mag to `galsim`'s"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"deltazp = 2.5 * np.log10(collecting_area)"
]
},
{
"cell_type": "markdown",
"metadata": {},
Expand All @@ -392,13 +420,13 @@
"source": [
"def sersic_mod(cutout, n, hlr, flux, pa, x0, y0, q):\n",
" nx, ny = cutout.array.shape\n",
" flux = flux*zp # total flux in photons/s\n",
" flux = flux*zp*deltazp # total flux in photons/s/cm^2\n",
" ser = galsim.Sersic(n, half_light_radius=hlr, flux=flux)\n",
" ser.shear(q=q, beta=pa*galsim.degrees) # change PA and axis ratio\n",
" offset = galsim.PositionD(x=x0, y=y0) # potentially shift a bit the profile\n",
" ser = galsim.convolve.Convolve(ser, psf_obj) # add PSF\n",
" img = galsim.ImageD(nx, ny) \n",
" ser.drawImage(image=img, offset=offset) # draw image of the model\n",
" ser = galsim.convolve.Convolve(ser, psf_obj) # add PSF -- the PSF is super-sampled!\n",
" img = galsim.ImageD(nx, ny, scale=0.11) \n",
" ser.drawImage(image=img, offset=offset, scale=0.11) # draw image of the model\n",
" return img.array.flatten() # return a 1D-array"
]
},
Expand All @@ -415,7 +443,7 @@
"metadata": {},
"outputs": [],
"source": [
"p0 = [catalog['n'][gal_sel][igal], catalog['half_light_radius'][gal_sel][igal], catalog['F158'][gal_sel][igal],\n",
"p0 = [catalog['n'][gal_sel][igal], catalog['half_light_radius'][gal_sel][igal], catalog['F106'][gal_sel][igal],\n",
" catalog['pa'][gal_sel][igal], 0., 0., catalog['ba'][gal_sel][igal]]\n",
"fid_model = sersic_mod(cutout, p0[0], p0[1], p0[2], p0[3], p0[4], p0[5], p0[6])"
]
Expand All @@ -433,7 +461,7 @@
"metadata": {},
"outputs": [],
"source": [
"err_img = img_asdf['roman']['err'][ymin:ymax, xmin:xmax].value.flatten()"
"err_img = img_asdf['err'][ymin:ymax, xmin:xmax].value.flatten()"
]
},
{
Expand Down Expand Up @@ -484,7 +512,7 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's see the differences"
"Let's examine the fit's residuals:"
]
},
{
Expand All @@ -494,7 +522,8 @@
"outputs": [],
"source": [
"plt.imshow(cutout.array-model_out.reshape(cutout.array.shape), origin='lower')\n",
"plt.colorbar()\n"
"plt.colorbar(label='DN/s')\n",
"plt.title('Fit residual');"
]
},
{
Expand All @@ -510,7 +539,7 @@
"metadata": {},
"outputs": [],
"source": [
"plt.imshow(fid_model.reshape(cutout.array.shape))\n",
"plt.imshow(fid_model.reshape(cutout.array.shape), norm=LogNorm(1, 3))\n",
"plt.colorbar(label='DN/s')\n",
"plt.title('Fiducial model')\n",
"plt.figure()\n",
Expand Down Expand Up @@ -545,8 +574,8 @@
"source": [
"## About this notebook\n",
"\n",
"**Author:** Javier Sánchez, Amethyst Barnes, Ami Choi. \n",
"**Updated On:** 2024-05-06"
"**Author:** Javier Sánchez ([email protected]), Amethyst Barnes, Ami Choi. \n",
"**Updated On:** 2024-06-17"
]
},
{
Expand All @@ -567,9 +596,9 @@
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"display_name": "Roman Calibration latest (2024-03-25)",
"language": "python",
"name": "python3"
"name": "roman-cal"
},
"language_info": {
"codemirror_mode": {
Expand Down
Loading
Loading