diff --git a/notebooks/romanisim_romancal/requirements.txt b/notebooks/romanisim_romancal/requirements.txt new file mode 100644 index 0000000..2a4573f --- /dev/null +++ b/notebooks/romanisim_romancal/requirements.txt @@ -0,0 +1,5 @@ +numpy +astropy +asdf +romancal +romanisim diff --git a/notebooks/romanisim_romancal/romanisim_romancal.ipynb b/notebooks/romanisim_romancal/romanisim_romancal.ipynb index b466c81..f28e283 100644 --- a/notebooks/romanisim_romancal/romanisim_romancal.ipynb +++ b/notebooks/romanisim_romancal/romanisim_romancal.ipynb @@ -32,13 +32,17 @@ "source": [ "## Imports\n", " Libraries used\n", + "- *subprocess* for running scripts from command line \n", "- *numpy* to handle array functions\n", "- *pandas* for working with tabular data\n", "- *matplotlib.pyplot* for plotting data\n", "- *astropy.io* wrting csv files\n", + "- *astropy.time* creating a time object \n", "- *astroquery.gaia* querying and downloading the Gaia data\n", - "- *romancal* for running the processing pipeline\n", - "- *asdf* for reading and writing asdf files" + "- *romanisim.gaia* querying and downloading the Gaia data\n", + "- *romancal.pipeline* for running the processing pipeline\n", + "- *asdf* for reading and writing asdf files\n", + "- *roman_datamodels* for reading and writing asdf files" ] }, { @@ -51,18 +55,19 @@ }, "outputs": [], "source": [ - "%matplotlib ipympl\n", - "#%matplotlib inline\n", + "%matplotlib inline\n", "import os\n", - "import sys\n", + "import subprocess\n", "import numpy as np\n", - "import pandas\n", "import matplotlib.pyplot as plt\n", - "from matplotlib.colors import LogNorm\n", + "from matplotlib.colors import LogNorm\n", "import astropy.io \n", "from astroquery.gaia import Gaia\n", "import asdf\n", - "from romancal.pipeline import ExposurePipeline" + "import roman_datamodels as rdm\n", + "from romancal.pipeline import ExposurePipeline\n", + "import romanisim.gaia\n", + "import astropy.time " ] }, { @@ -129,11 +134,12 @@ }, "outputs": [], "source": [ - "def make_catalog_synthetic(filename, ra, dec, bandpass='F158', n_ser=0, n_psf=0, area=0.24, seed=12):\n", + "\n", + "def make_catalog_synthetic(outfile_name, ra, dec, bandpass='F158', n_ser=0, n_psf=0, area=0.24, seed=12):\n", " \"\"\"\n", " Input\n", " -----\n", - " fileanme (str): output file name\n", + " outfile_name (str): output file name\n", " ra (float): ra of the center of WFI [deg]\n", " dec (float): dec of the center of WFI [deg]\n", " bandpass (str): Photometric bandpasss e.g F158\n", @@ -143,38 +149,39 @@ " seed (int): Seed for random number generator\n", " \"\"\"\n", " np.random.seed(seed)\n", - " scale=np.sqrt(area/18.0)\n", + " scale = np.sqrt(area/18.0)\n", " # total number of sources\n", - " n_sources=n_psf+n_ser\n", - " d={}\n", - " d['ra']=ra+(np.random.uniform(size=n_sources)-0.5)*scale*6/np.cos(np.radians(dec))\n", - " d['dec']=dec+(np.random.uniform(size=n_sources)-0.5)*scale*3 \n", - " # set the type 'PSF' and 'SER' \n", - " d['type']=np.array(['PSF']*n_psf+['SER']*n_ser)\n", + " n_sources = n_psf+n_ser\n", + " d = {}\n", + " d['ra'] = ra+(np.random.uniform(size=n_sources) - 0.5) * \\\n", + " scale*6/np.cos(np.radians(dec))\n", + " d['dec'] = dec+(np.random.uniform(size=n_sources)-0.5)*scale*3\n", + " # set the type 'PSF' and 'SER'\n", + " d['type'] = np.array(['PSF']*n_psf+['SER']*n_ser)\n", " # set sersic index of sources in range [2, 4]\n", - " d['n']=np.int64(2+np.random.uniform(size=n_sources)*2)\n", + " d['n'] = np.int64(2+np.random.uniform(size=n_sources)*2)\n", " # set half light radius radius [3, 8]\n", - " d['half_light_radius']=3+np.random.uniform(size=n_sources)*5\n", + " d['half_light_radius'] = 3+np.random.uniform(size=n_sources)*5\n", " # set position angle in range [0,90]\n", - " d['pa']=np.random.uniform(size=n_sources)*90.0\n", + " d['pa'] = np.random.uniform(size=n_sources)*90.0\n", " # set major to minor axis ratio in range [0.3,0.7]\n", - " d['ba']=0.3+np.random.uniform(size=n_sources)*0.7\n", - " # set bandpass flux \n", - " temp1=np.power(10.0,np.random.uniform(size=n_psf))*1e-8\n", - " temp2=np.power(10.0,np.random.uniform(size=n_ser))*1e-7\n", - " d['F158']=np.concatenate([temp1,temp2])\n", + " d['ba'] = 0.3+np.random.uniform(size=n_sources)*0.7\n", + " # set bandpass flux\n", + " temp1 = np.power(10.0, np.random.uniform(size=n_psf))*1e-8\n", + " temp2 = np.power(10.0, np.random.uniform(size=n_ser))*1e-7\n", + " d['F158'] = np.concatenate([temp1, temp2])\n", " # write the catalog file\n", - " #pandas.DataFrame(d).to_csv(filename,sep=' ',index=False)\n", - " astropy.io.ascii.write(d,filename,overwrite=True,names=['ra','dec','type','n','half_light_radius','pa','ba','F158']) \n", - " \n", + " astropy.io.ascii.write(d, outfile_name, overwrite=True, names=[\n", + " 'ra', 'dec', 'type', 'n', 'half_light_radius', 'pa', 'ba', 'F158'])\n", "\n", - "dirname='test_data/'\n", - "if not os.path.exists(dirname):\n", - " os.makedirs(dirname)\n", "\n", + "# dirname = 'test_data/'\n", + "# if not os.path.exists(dirname):\n", + "# os.makedirs(dirname)\n", + "dirname = ''\n", "ra_cen, dec_cen = (0.0, 0.0)\n", - "bandpass='F158'\n", - "make_catalog_synthetic(dirname+'catalog_mixed.csv', ra_cen, dec_cen, n_ser=900,n_psf=360)\n" + "make_catalog_synthetic(dirname+'catalog_mixed.csv',\n", + " ra_cen, dec_cen, n_ser=900, n_psf=360)" ] }, { @@ -182,7 +189,7 @@ "metadata": {}, "source": [ "### Make input catalog for romanisim from Gaia\n", - "We query the Gaia archive to do a cone serach around user specified equatorial coordiantes and synthetic photometry from table ```gaiadr3.synthetic_photometry_gspc``` is added to it. Flux in roman bandpass is estimated assuming that it is 1/100 th of the flux in filter ```f814w_acswfc_mag```. " + "We query the Gaia archive to do a cone serach around user specified equatorial coordiantes and then wriyte it to an [```.ecsv```](https://docs.astropy.org/en/stable/io/ascii/ecsv.html) file (a csv file with extra commented out rows at the top providing details about the columns). Flux in roman bandpass is assumed to be 1/100 times (5 magnitudes) fainter than flux in gaia bandpass ```phot_g_mean_mag```." ] }, { @@ -192,35 +199,31 @@ "outputs": [], "source": [ "\n", - "def make_catalog_gaia(filename, ra, dec, bandpass='F158', radius=1.0, row_limit=1000):\n", - " gaia_filter_name='f814w_acswfc_mag'\n", + "def make_catalog_gaia(outfile_name, ra, dec, date='2026-01-01', bandpasses=['F158'], radius=1.0, row_limit=1000):\n", + " \"\"\"\n", + " date (str): 'YYYY-MM-DD' e.g. '2027-06-16' \n", + " radius (float): search radius in degrees\n", + " row_limit (int): Use -1 for unlimited\n", + " bandpasses list(str): list of roman bandpasses. If set to None the full list is loaded as specified in. \n", + " \"\"\"\n", " row_limit = f'TOP {row_limit}' if row_limit > -1 else ''\n", - " query = f\"SELECT {row_limit} source.source_id, source.ra, source.ra_error, \" \\\n", - " f\"source.dec, source.dec_error, synphot.{gaia_filter_name} \" \\\n", - " f\"FROM gaiadr3.gaia_source AS source, \" \\\n", - " f\"gaiadr3.synthetic_photometry_gspc AS synphot \" \\\n", - " f\"WHERE 1 = CONTAINS(POINT('ICRS', source.ra, source.dec), \" \\\n", - " f\"CIRCLE('ICRS',{ra}, {dec}, {radius})) \" \\\n", - " f\"AND source.source_id = synphot.source_id\"\n", + " q = f'select {row_limit} * from gaiadr3.gaia_source where distance({ra}, {dec}, ra, dec) < {radius}'\n", + " job = Gaia.launch_job_async(q)\n", + " r = job.get_results()\n", + "\n", + " r['phot_g_mean_mag'] = r['phot_g_mean_mag']+5\n", + " if bandpasses is None:\n", + " bandpasses = set(romanisim.bandpass.galsim2roman_bandpass.values())\n", "\n", - " job = Gaia.launch_job_async(query)\n", - " print(job)\n", - " results = job.get_results()\n", + " if outfile_name.endswith('.ecsv'):\n", + " romanisim.gaia.gaia2romanisimcat(r, astropy.time.Time(\n", + " date), fluxfields=set(bandpasses)).write(outfile_name, overwrite=True)\n", + " else:\n", + " raise RuntimeError('output file name should end with .ecsv')\n", "\n", - " d={}\n", - " n_sources=results['ra'].size\n", - " d['ra']=results['ra']\n", - " d['dec']=results['dec']\n", - " d['type']=np.array(['PSF']*n_sources)\n", - " d['n']=np.array([2.0]*n_sources)\n", - " d['half_light_radius']=np.array([2.0]*n_sources)\n", - " d['pa']=np.array([90.0]*n_sources)\n", - " d['ba']=np.array([0.3]*n_sources)\n", - " d['F158']=np.power(10.0,-results['f814w_acswfc_mag']/2.5)/100\n", - " \n", - " astropy.io.ascii.write(d,filename,overwrite=True,names=['ra','dec','type','n','half_light_radius','pa','ba','F158']) \n", "\n", - "make_catalog_gaia(f\"{dirname}catalog_gaia.csv\",ra_cen, dec_cen)" + "# dirname='test_data/'\n", + "make_catalog_gaia(f\"{dirname}catalog_gaia.ecsv\", ra_cen, dec_cen)" ] }, { @@ -237,12 +240,8 @@ "metadata": {}, "outputs": [], "source": [ - "import subprocess\n", - "print(subprocess.getoutput('romanisim-make-image -h'))\n", - "#result=subprocess.run(['romanisim-make-image', '-h'], capture_output=True,text=True)\n", - "#print(result.stdout)\n", - "# This suppresses the stdout\n", - "#os.system('romanisim-make-image -h')" + "\n", + "print(subprocess.getoutput('romanisim-make-image -h'))" ] }, { @@ -258,24 +257,24 @@ "metadata": {}, "outputs": [], "source": [ - "sca=2\n", - "exposure_no=1\n", - "date='2026-01-01T01:00:00'\n", - "dirname='test_data/'\n", - "l1_fname=f\"romanisim_mixed_sca{sca}\"\n", - "catalog_filename=dirname+'catalog_mixed.csv'\n", + "sca = 2\n", + "exposure_no = 1\n", + "date = '2026-01-01T01:00:00'\n", + "dirname = 'test_data/'\n", + "l1_fname = f\"romanisim_mixed_sca{sca}\"\n", + "catalog_filename = dirname+'catalog_mixed.csv'\n", "\n", - "level=1\n", - "filename=f\"{dirname}{l1_fname}_uncal.asdf\" \n", + "level = 1\n", + "filename = f\"{dirname}{l1_fname}_uncal.asdf\"\n", "cmd = f\"romanisim-make-image --radec {ra_cen} {dec_cen} --catalog {catalog_filename} --usecrds --webbpsf --level {level} --bandpass F158 --date {date} --rng_seed 11 --sca {sca} {filename}\"\n", "print(cmd)\n", "os.system(cmd)\n", "\n", - "level=2\n", - "filename=f\"{dirname}{l1_fname}_l2.asdf\" \n", + "level = 2\n", + "filename = f\"{dirname}{l1_fname}_l2.asdf\"\n", "cmd = f\"romanisim-make-image --radec {ra_cen} {dec_cen} --catalog {catalog_filename} --usecrds --webbpsf --level {level} --bandpass F158 --date {date} --rng_seed 11 --sca {sca} {filename}\"\n", "print(cmd)\n", - "os.system(cmd)\n" + "os.system(cmd)" ] }, { @@ -292,17 +291,18 @@ "metadata": {}, "outputs": [], "source": [ - "l1_fname='romanisim_mixed_sca2'\n", - "filename=f\"{dirname}{l1_fname}_uncal.asdf\" \n", - "result = ExposurePipeline.call(filename,save_results=True,output_dir=dirname,\n", + "l1_fname = 'romanisim_mixed_sca2'\n", + "filename = f\"{dirname}{l1_fname}_uncal.asdf\"\n", + "result = ExposurePipeline.call(filename, save_results=True, output_dir=dirname,\n", " steps={'source_detection':\n", - " {'save_catalogs':True},\n", - " 'rampfit':{'threshold_intercept':5.5,\n", - " 'threshold_constant':0.33}})\n", + " {'save_catalogs': True},\n", + " 'rampfit': {'threshold_intercept': 5.5,\n", + " 'threshold_constant': 0.33}})\n", "\n", "# move catalog file to data directory\n", "if os.path.isfile(f\"{l1_fname}_uncal_tweakreg_catalog.asdf\"):\n", - " os.system(f'mv {l1_fname}_uncal_tweakreg_catalog.asdf {dirname}{l1_fname}_uncal_tweakreg_catalog.asdf')\n" + " os.system(\n", + " f'mv {l1_fname}_uncal_tweakreg_catalog.asdf {dirname}{l1_fname}_uncal_tweakreg_catalog.asdf')" ] }, { @@ -310,7 +310,9 @@ "metadata": {}, "source": [ "## Open and explore the L1 and L2 images generated by romanisim \n", - "We open the Level 1 file, its cooresponding Level 2 file and the catalog sources identfied in the Level 2 file. " + "We open the Level 1 file, its cooresponding Level 2 file and the catalog of sources identfied in the Level 2 file. We use ```roman_datamodels``` to open the first two files as it provides attribute style access to the underlying data, but one can also open them ```AsdfObject``` object for dictionary style access (see notebook on How to open Roman Data Files for more details). \n", + "\n", + "We examine the shape and the units of the data in the files. The L1 files are of shape 6x4096x4096 and have units of DN. The first dimension represents the time. The L2 files are of shape 4088x4088 and are in units of electrons/s. There is no time dimension as we have computed the slope or rate. The shape of image is smaller by 8 pixels as the 4 pixel border which has no science data has been removed. Finally, we print the sources contained in the catalog. The x and y coordinates of the centroid and the flux are shown. The centroid is in units of pixel and ranges between 0 to 4088." ] }, { @@ -319,35 +321,40 @@ "metadata": {}, "outputs": [], "source": [ + "\n", "# read uncalibrated data\n", - "l1_fname='romanisim_mixed_sca2'\n", - "af1=asdf.open(f\"{dirname}{l1_fname}_uncal.asdf\",copy_arrays=True)\n", + "l1_fname = 'romanisim_mixed_sca2'\n", + "rdm1 = rdm.open(f\"{dirname}{l1_fname}_uncal.asdf\")\n", "\n", "# read calibrated data\n", - "af2=asdf.open(f\"{dirname}{l1_fname}_cal.asdf\",copy_arrays=True)\n", + "rdm2 = rdm.open(f\"{dirname}{l1_fname}_cal.asdf\")\n", "\n", "# Read the source catalog\n", - "catalog = asdf.open(f\"{dirname}{l1_fname}_uncal_tweakreg_catalog.asdf\")['tweakreg_catalog'].copy()\n", + "catalog = asdf.open(f\"{dirname}{l1_fname}_uncal_tweakreg_catalog.asdf\")[\n", + " 'tweakreg_catalog'].copy()\n", "\n", "\n", "# explore shape and units of l1 and l2 data\n", - "print(f\"l1 data shape: {str(af1['roman']['data'].shape):16} units: {af1['roman']['data'].unit}\")\n", - "print(f\"l2 data shape: {str(af2['roman']['data'].shape):16} units: {af2['roman']['data'].unit}\")\n", + "print(\n", + " f\"l1 data shape: {str(af1['roman']['data'].shape):16} units: {af1['roman']['data'].unit}\")\n", + "print(\n", + " f\"l2 data shape: {str(af2['roman']['data'].shape):16} units: {af2['roman']['data'].unit}\")\n", "print()\n", "\n", "# print the source catalog\n", - "catalog = asdf.open(f\"{dirname}{l1_fname}_uncal_tweakreg_catalog.asdf\")['tweakreg_catalog'].copy()\n", "print(f\"{'index':6} {'xcentroid':8} {'ycentroid':8} {'flux':8}\")\n", "print(\"-------------------------------------\")\n", "for i in range(catalog.size):\n", - " print(f\"{i:<6} {catalog['xcentroid'][i]:<8.2f} {catalog['ycentroid'][i]:8.2f} {catalog['flux'][i]:<8.2f}\")\n" + " print(\n", + " f\"{i:<6} {catalog['xcentroid'][i]:<8.2f} {catalog['ycentroid'][i]:8.2f} {catalog['flux'][i]:<8.2f}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "### Print information contained in the files" + "### Print information contained in the files\n", + "The information in the roman_datamodel object is stored in a hierarchical fashion and can be accessed using the attribute notation. The L1 object has three high level keys ```meta, data``` and ```amp33```, the L2 object has many more. To view the full list of available attributes and subattributes we use ```.info()``` member function. Try to locate the following atrribiutes, which we will use later, ```rdm1.meta.exposure.frame_time``` and ```rdm1.meta.exposure.read_pattern``` . " ] }, { @@ -356,8 +363,21 @@ "metadata": {}, "outputs": [], "source": [ - "af1.info(max_rows=None)\n", - "af2.info(max_rows=None)" + "print('rdm1 keys')\n", + "print('-----------')\n", + "for key in rdm1.keys():\n", + " print(key)\n", + "print('-----------')\n", + "print('rdm2 keys')\n", + "for key in rdm2.keys():\n", + " print(key)\n", + "print('-----------')\n", + "print('rdm1 info')\n", + "rdm1.info(max_rows=None)\n", + "print('-----------')\n", + "print('rdm2 info')\n", + "rdm2.info(max_rows=None)\n", + "print('-----------')" ] }, { @@ -365,7 +385,8 @@ "metadata": {}, "source": [ "### View the Level 2 image \n", - "All panels correspond to the same image but their extent are different. Each panel is a zoomed in version of the previous one." + "We now view the Level 2 images. All panels correspond to the same image but their extent are different. Each panel is a zoomed in version of the previous one.\n", + "One can see some points source objects (stars descibed by a PSF) and extended source objects (galaxies described by a sersic profile). The white spots are bad pixels and are to be ignored. The objects can be seen more clearly in the zoomed bottom left panel. One can read the coordinates and the flux values by sliding the cursor on the panels.\n" ] }, { @@ -375,25 +396,26 @@ "outputs": [], "source": [ "\n", - "image=af2['roman']['data'].value\n", - "fig, axs = plt.subplots(2,2)\n", - "i,j=(4088//2,4088//2)\n", - "norm = LogNorm(vmin=2,vmax=20)\n", - "im1=axs.flat[0].imshow(image,cmap='YlOrBr_r',norm=norm,origin='lower')\n", - "im1=axs.flat[1].imshow(image,cmap='YlOrBr_r',norm=norm,origin='lower')\n", - "axs.flat[1].axis([i-1024,i+1024,j-1024,j+1024])\n", - "im1=axs.flat[2].imshow(image,cmap='YlOrBr_r',norm=norm,origin='lower')\n", - "axs.flat[2].axis([i-256,i+256,j-256,j+256])\n", - "im1=axs.flat[3].imshow(image,cmap='YlOrBr_r',norm=norm,origin='lower')\n", - "axs.flat[3].axis([i-64,i+64,j-64,j+64])\n", - "fig.colorbar(im1,ax=axs)\n" + "image = rdm2.data.value\n", + "fig, axs = plt.subplots(2, 2)\n", + "i, j = (4088//2, 4088//2)\n", + "norm = LogNorm(vmin=2, vmax=20)\n", + "im1 = axs.flat[0].imshow(image, cmap='YlOrBr_r', norm=norm, origin='lower')\n", + "im1 = axs.flat[1].imshow(image, cmap='YlOrBr_r', norm=norm, origin='lower')\n", + "axs.flat[1].axis([i-1024, i+1024, j-1024, j+1024])\n", + "im1 = axs.flat[2].imshow(image, cmap='YlOrBr_r', norm=norm, origin='lower')\n", + "axs.flat[2].axis([i-256, i+256, j-256, j+256])\n", + "im1 = axs.flat[3].imshow(image, cmap='YlOrBr_r', norm=norm, origin='lower')\n", + "axs.flat[3].axis([i-64, i+64, j-64, j+64])\n", + "fig.colorbar(im1, ax=axs)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "### View 64x64 pixel cutouts of Level 2 images centered on detected sources" + "### View 64x64 pixel cutouts of Level 2 images centered on detected sources\n", + "We choose four sources from the source catalog and plot 64x64 pixel image cutouts centered on them. The first and fourth are point sources, while the second and third are extended sources. One can see that the point sources are tightly confined due to the PSF being very narrow and is undersampled. " ] }, { @@ -403,14 +425,17 @@ "outputs": [], "source": [ "# plot 64x64 cutouts of l2 images centered on detected sources\n", - "norm=LogNorm(vmin=2,vmax=100)\n", - "fig, axs = plt.subplots(2,2)\n", - "for j,i in enumerate([5, 6, 7 , 8]):\n", - " width=32\n", - " jp=int(catalog['xcentroid'][i])\n", - " ip=int(catalog['ycentroid'][i]) \n", - " im1=axs.flat[j].imshow(image[ip-width:ip+width,jp-width:jp+width],cmap='YlOrBr_r',norm=norm,origin='lower') \n", - "fig.colorbar(im1,ax=axs)\n" + "norm = LogNorm(vmin=2, vmax=100)\n", + "fig, axs = plt.subplots(2, 2)\n", + "for j, i in enumerate([5, 6, 7, 8]):\n", + " width = 32\n", + " jp = int(catalog['xcentroid'][i])\n", + " ip = int(catalog['ycentroid'][i])\n", + " im1 = axs.flat[j].imshow(image[ip-width:ip+width, jp-width:jp+width],\n", + " cmap='YlOrBr_r', norm=norm, origin='lower')\n", + " axs.flat[j].set_title('Source No '+str(i))\n", + "fig.tight_layout()\n", + "fig.colorbar(im1, ax=axs)" ] }, { @@ -418,10 +443,9 @@ "metadata": {}, "source": [ "### Plot ramps centered on detected sources from Level 1 data and compare its slope with the value in the Level 2 data \n", - "To plot the ramps, we being by first computing the mean time of each resultant from the read out pattern. We \n", - "do crude estimate of slope by taking the difference between the first and the last resultant values and dividing it by the time difference.\n", + "To plot the ramps, we being by first computing the mean time of each resultant from the read out pattern. We do crude estimate of slope by taking the difference between the first and the last resultant values and dividing it by the time difference.\n", "\n", - "A few things should be kept in mind while working with Level 1 and Level 2 data. The shape of Level 1 and Level 2 are different, as there exists a 4 pixel wide border which is ignored in Level 2. Hence, a pixel (i,j) in Level 2 corresponds to pixel (i+4,j+4) in Level 1. Also the units of Level 1 are in DN while that of Level 2 data are electron/s. The ratio of electrons to DN is know as the gain which can be found in the calibration reference file." + "A few things should be kept in mind while working with Level 1 and Level 2 data. The shape of Level 1 and Level 2 are different, as there exists a 4 pixel wide border which is ignored in Level 2. Hence, a pixel (i,j) in Level 2 corresponds to pixel (i+4,j+4) in Level 1. Also the units of Level 1 are in DN while that of Level 2 data are electron/s. The ratio of electrons to DN is know as the gain which can be found in the calibration reference file. From the printed results we see that our estimated slope is in good agreement with what is shown in Level 2 data except for a factor of 2, which is due to the gain. " ] }, { @@ -431,30 +455,33 @@ "outputs": [], "source": [ "\n", + "frame_time = rdm1.meta.exposure.frame_time\n", + "read_pattern = rdm1.meta.exposure.read_pattern\n", "# mean time of each resultant\n", - "frame_time=af1['roman']['meta']['exposure']['frame_time']\n", - "read_pattern=af1['roman']['meta']['exposure']['read_pattern']\n", - "t_mean=np.array([frame_time*np.mean(t) for t in read_pattern])\n", - "image1=af1['roman']['data'].value\n", - "image2=af2['roman']['data'].value\n", + "t_mean = np.array([frame_time*np.mean(t) for t in read_pattern])\n", + "image1 = rdm1.data.value\n", + "image2 = rdm2.data.value\n", "\n", - "fig,ax=plt.subplots()\n", - "print(f\"{'source no':16} {'Level 2 [i,j]':16} {'ramp slope':16} {'Level 2 data':16}\")\n", - "print(f\"{'':16} {'':16} {str(af1['roman']['data'].unit)+'/s':16} {af2['roman']['data'].unit:16}\")\n", + "fig, ax = plt.subplots()\n", + "print(\n", + " f\"{'source no':16} {'Level 2 [i,j]':16} {'ramp slope':16} {'Level 2 data':16}\")\n", + "print(\n", + " f\"{'':16} {'':16} {str(af1['roman']['data'].unit)+'/s':16} {af2['roman']['data'].unit:16}\")\n", "print(\"---------------------------------------------------------------------\")\n", "for k in [3, 4, 5, 6]:\n", " # image coordinates in Level 1\n", - " ip,jp=(int(catalog['ycentroid'][k])+4,int(catalog['xcentroid'][k])+4)\n", - " ax.plot(t_mean,image1[:,ip,jp],'-',label=f\"source no: {k}, L2_image[{ip},{jp}]\")\n", - " ax.plot(t_mean,image1[:,ip,jp],'k.')\n", + " ip, jp = (int(catalog['ycentroid'][k])+4, int(catalog['xcentroid'][k])+4)\n", + " ax.plot(t_mean, image1[:, ip, jp], '-',\n", + " label=f\"source no: {k}, L2_image[{ip},{jp}]\")\n", + " ax.plot(t_mean, image1[:, ip, jp], 'k.')\n", " # crude slope from Level 1\n", - " dt=t_mean[-1]-t_mean[0]\n", - " slope=(image1[-1,ip,jp]-image1[0,ip,jp])/dt\n", - " print(f\"{k:<16} [{ip:<4},{jp:<4}]{'':5} {slope:<16.2f} {image2[ip-4,jp-4]:<16.2f}\")\n", + " dt = t_mean[-1]-t_mean[0]\n", + " slope = (image1[-1, ip, jp]-image1[0, ip, jp])/dt\n", + " print(\n", + " f\"{k:<16} [{ip:<4}, {jp:<4}]{'':5} {slope:<16.2f} {image2[ip-4, jp-4]:<16.2f}\")\n", "ax.set_xlabel('Mean resultant time [s]')\n", "ax.set_ylabel(f\"counts [{af1['roman']['data'].unit}]\")\n", - "ax.legend()\n", - "\n" + "ax.legend()" ] }, { @@ -477,7 +504,7 @@ "source": [ "## About this notebook\n", "**Author:** Sanjib Sharma, Roman Telescope Branch. \n", - "**Updated On:** 2024-05-01" + "**Updated On:** 2024-23-01" ] }, { @@ -519,7 +546,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.11.9" + "version": "3.9.12" } }, "nbformat": 4,