Skip to content

Commit

Permalink
Merge pull request #31 from fjaviersanchez/js_s3updates
Browse files Browse the repository at this point in the history 
updated with s3df file stream
  • Loading branch information
@tddesjardins
tddesjardins authored Jun 17, 2024
2 parents 120b2f6 + 8f1bdec commit 2f34e82
Show file tree
Hide file tree
Showing 2 changed files with 300 additions and 126 deletions.
113 changes: 71 additions & 42 deletions content/notebooks/measuring_galaxy_shapes/measuring_galaxy_shapes.ipynb
View file
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

0 comments on commit 2f34e82

Please sign in to comment.